From f557d1f49a2914c6084dd18efc783395228d8ce0 Mon Sep 17 00:00:00 2001 From: Paul Duncan Date: Tue, 5 Feb 2019 00:22:15 -0500 Subject: mv *.[hc] src/ --- .gitignore | 2 +- Makefile | 11 +- km-draw.c | 31 - km-find.c | 115 ---- km-init-forgy.c | 51 -- km-init-kmeans.c | 124 ---- km-init-rand.c | 43 -- km-init.c | 85 --- km-load.c | 112 ---- km-rand-erand48.c | 89 --- km-rand-libc.c | 56 -- km-rand-path.c | 152 ----- km-rand.c | 33 -- km-set.c | 267 --------- km-solve.c | 200 ------- km.h | 230 -------- main.c | 398 ------------- src/km-draw.c | 31 + src/km-find.c | 115 ++++ src/km-init-forgy.c | 51 ++ src/km-init-kmeans.c | 124 ++++ src/km-init-rand.c | 43 ++ src/km-init.c | 85 +++ src/km-load.c | 112 ++++ src/km-rand-erand48.c | 89 +++ src/km-rand-libc.c | 56 ++ src/km-rand-path.c | 152 +++++ src/km-rand.c | 33 ++ src/km-set.c | 267 +++++++++ src/km-solve.c | 200 +++++++ src/km.h | 230 ++++++++ src/main.c | 398 +++++++++++++ src/stb_image_write.h | 1568 +++++++++++++++++++++++++++++++++++++++++++++++++ src/util.h | 49 ++ stb_image_write.h | 1568 ------------------------------------------------- util.h | 49 -- 36 files changed, 3610 insertions(+), 3609 deletions(-) delete mode 100644 km-draw.c delete mode 100644 km-find.c delete mode 100644 km-init-forgy.c delete mode 100644 km-init-kmeans.c delete mode 100644 km-init-rand.c delete mode 100644 km-init.c delete mode 100644 km-load.c delete mode 100644 km-rand-erand48.c delete mode 100644 km-rand-libc.c delete mode 100644 km-rand-path.c delete mode 100644 km-rand.c delete mode 100644 km-set.c delete mode 100644 km-solve.c delete mode 100644 km.h delete mode 100644 main.c create mode 100644 src/km-draw.c create mode 100644 src/km-find.c create mode 100644 src/km-init-forgy.c create mode 100644 src/km-init-kmeans.c create mode 100644 src/km-init-rand.c create mode 100644 src/km-init.c create mode 100644 src/km-load.c create mode 100644 src/km-rand-erand48.c create mode 100644 src/km-rand-libc.c create mode 100644 src/km-rand-path.c create mode 100644 src/km-rand.c create mode 100644 src/km-set.c create mode 100644 src/km-solve.c create mode 100644 src/km.h create mode 100644 src/main.c create mode 100644 src/stb_image_write.h create mode 100644 src/util.h delete mode 100644 stb_image_write.h delete mode 100644 util.h diff --git a/.gitignore b/.gitignore index fa86ed6..b54b14d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,3 @@ -*.o +src/*.o km-test *.png diff --git a/Makefile b/Makefile index 0372ecf..19691c7 100644 --- a/Makefile +++ b/Makefile @@ -1,8 +1,9 @@ APP=km-test CFLAGS=-W -Wall -Wextra -Werror -pedantic -std=c11 -O2 -OBJS=km-set.o km-draw.o km-load.o km-find.o km-rand.o km-solve.o \ - km-init-rand.o km-init-forgy.o km-init-kmeans.o km-init.o main.o \ - km-rand-libc.o km-rand-path.o km-rand-erand48.o +OBJS=src/km-set.o src/km-draw.o src/km-load.o src/km-find.o \ + src/km-rand.o src/km-solve.o src/km-init-rand.o src/main.o \ + src/km-init-forgy.o src/km-init-kmeans.o src/km-init.o \ + src/km-rand-libc.o src/km-rand-path.o src/km-rand-erand48.o LIBS=-lm .PHONY=all clean @@ -11,8 +12,8 @@ all: $(APP) $(APP): $(OBJS) $(CC) $(CFLAGS) -o $(APP) $(OBJS) $(LIBS) -%.o: %.c +%.o: src/%.c $(CC) -c $(CFLAGS) $< clean: - $(RM) $(OBJS) + $(RM) $(APP) $(OBJS) diff --git a/km-draw.c b/km-draw.c deleted file mode 100644 index ddf1b26..0000000 --- a/km-draw.c +++ /dev/null @@ -1,31 +0,0 @@ -#include // size_t -#include "util.h" -#include "km.h" - -void -km_set_draw( - const km_set_t * const set, - uint8_t * const rgb, - const size_t width, - const size_t height, - const int dot_size, - const uint32_t color -) { - for (size_t i = 0; i < set->num_rows; i++) { - const float *row = km_set_get_row(set, i); - - for (int yo = 0; yo < dot_size; yo++) { - for (int xo = 0; xo < dot_size; xo++) { - const int x = (width - 1) * row[0] - dot_size / 2 + xo, - y = (height - 1) * row[1] - dot_size / 2 + yo, - ofs = 3 * (width * y + x); - - if (x >= 0 && x < (int) width && y >= 0 && y < (int) height) { - rgb[ofs + 0] = (color & 0xff0000) >> 16; - rgb[ofs + 1] = (color & 0x00ff00) >> 8; - rgb[ofs + 2] = (color & 0x0000ff); - } - } - } - } -} diff --git a/km-find.c b/km-find.c deleted file mode 100644 index 657c645..0000000 --- a/km-find.c +++ /dev/null @@ -1,115 +0,0 @@ -#include // bool -#include // fabsf() -#include "util.h" -#include "km.h" - -#define MIN_CLUSTER_DISTANCE 0.0001 - -static float -get_mean_cluster_size( - const km_set_t * const set -) { - float sum = 0; - size_t num_filled = 0; - - for (size_t i = 0; i < set->num_rows; i++) { - if (set->ints[i] > 0) { - sum += set->ints[i]; - num_filled++; - } - } - - return (num_filled > 0) ? (sum / num_filled) : 0; -} - -static void -find_solve_on_stats( - const km_set_t * const set, - const km_solve_stats_t * const stats, - void * const cb_data -) { - km_find_data_t * const find_data = cb_data; - UNUSED(set); - - // save total sum and silouette - find_data->distance_sum = stats->sum; - find_data->silouette = stats->silouette; -} - -static const km_solve_cbs_t -FIND_SOLVE_CBS = { - .on_stats = find_solve_on_stats, -}; - -bool -km_find( - const km_set_t * const set, - const km_find_cbs_t * const cbs, - void * const cb_data -) { - // check init callback - if (!cbs->on_init) { - return false; - } - - // check fini callback - if (!cbs->on_fini) { - return false; - } - - // check data callback - if (!cbs->on_data) { - return false; - } - - float best_silouette = -2.0; - for (size_t i = 2; i < cbs->max_clusters; i++) { - for (size_t j = 0; j < cbs->num_tests; j++) { - // init cluster set - km_set_t cs; - if (!cbs->on_init(&cs, i, set, cb_data)) { - // return failure - return false; - } - - // init find data - km_find_data_t find_data = { - .cluster_set = &cs, - .num_clusters = i, - }; - - // solve test - // (populates sum and silouette) - if (!km_solve(&cs, set, &FIND_SOLVE_CBS, &find_data)) { - // return failure - return false; - } - - // populate mean cluster size - find_data.mean_cluster_size = get_mean_cluster_size(&cs); - - // emit result - cbs->on_data(&find_data, cb_data); - - if (find_data.silouette > best_silouette) { - // emit new best result - if (cbs->on_best && !cbs->on_best(find_data.silouette, &cs, cb_data)) { - // return failure - return false; - } - - // update best silouette - best_silouette = find_data.silouette; - } - - // finalize cluster set - if (!cbs->on_fini(&cs, cb_data)) { - // return failure - return false; - } - } - } - - // return success - return true; -} diff --git a/km-init-forgy.c b/km-init-forgy.c deleted file mode 100644 index 5ff0c12..0000000 --- a/km-init-forgy.c +++ /dev/null @@ -1,51 +0,0 @@ -#include // bool -#include // memset() -#include "util.h" -#include "km.h" - -// init a set with num_clusters clusters of shape num_floats by picking -// random initial points from the set -bool -km_init_forgy( - km_set_t * const cs, - const size_t num_clusters, - const km_set_t * const set, - km_rand_t * const rs -) { - const size_t num_floats = set->shape.num_floats, - stride = sizeof(float) * num_floats; - - // init cluster shape - const km_shape_t shape = { - .num_floats = num_floats, - .num_ints = 1, - }; - - // get random row offsets - size_t rows[num_clusters]; - if (!km_rand_get_sizes(rs, num_clusters, rows)) { - // return failure - return false; - } - - // generate random cluster centers - float floats[num_floats * num_clusters]; - for (size_t i = 0; i < num_clusters; i++) { - const size_t row_num = rows[i] % set->num_rows; - const float * const row_floats = km_set_get_row(set, row_num); - memcpy(floats + i * num_floats, row_floats, stride); - } - - // FIXME: should probably be heap-allocated - int ints[num_clusters]; - memset(ints, 0, sizeof(ints)); - - // init cluster set - if (!km_set_init(cs, &shape, num_clusters)) { - // return failure - return false; - } - - // add data, return result - return km_set_push(cs, num_clusters, floats, ints); -} diff --git a/km-init-kmeans.c b/km-init-kmeans.c deleted file mode 100644 index f0c1b5f..0000000 --- a/km-init-kmeans.c +++ /dev/null @@ -1,124 +0,0 @@ -#include // bool -#include // memset() -#include // FLT_MAX -#include "util.h" -#include "km.h" - -// sum the squared distance of every row in set from this point -static inline float -sum_distance_squared( - const float * const floats, - const km_set_t * const set -) { - float r = 0; - - // sum squared distances - for (size_t i = 0; i < set->num_rows; i++) { - const float * const vals = km_set_get_row(set, i); - r += distance_squared(set->shape.num_floats, floats, vals); - } - - // return result - return r; -} - -static const float * -get_random_row( - const km_set_t * const set, - km_rand_t * const rs -) { - // get random offset - size_t ofs = 0; - if (!km_rand_get_sizes(rs, 1, &ofs)) { - die("km_rand_fill_sizes()"); - } - - return km_set_get_row(set, ofs % set->num_rows); -} - -// init a cluster set of num_clusters by picking the first cluster from -// the set of points at random -// random initial points from the set -bool -km_init_kmeans( - km_set_t * const cs, - const size_t num_clusters, - const km_set_t * const set, - km_rand_t * const rs -) { - const size_t num_floats = set->shape.num_floats, - stride = sizeof(float) * num_floats; - - // row values (filled below) - float sums[num_clusters]; - float floats[num_floats * num_clusters]; - - // pick first row randomly - { - // get random row - const float * const vals = get_random_row(set, rs); - - // copy row values to floats buffer - memcpy(floats, vals, stride); - - // copy row values to floats buffer - sums[0] = sum_distance_squared(floats, set); - } - - for (size_t i = 1; i < num_clusters;) { - // get random row - const float * const vals = get_random_row(set, rs); - - // calculate squared distance to nearest cluster - size_t ofs = 0; - float min_d2 = FLT_MAX; - for (size_t j = 0; j < i; j++) { - const float d2 = distance_squared(num_floats, vals, floats + j * stride); - if (d2 < min_d2) { - ofs = j; - min_d2 = d2; - } - } - - // get a random floating point value - float rand_val = 0; - if (!km_rand_get_floats(rs, 1, &rand_val)) { - // return failure - return false; - } - - // check random value - if (rand_val / sums[ofs] < min_d2) { - // get destination - float * const dst = floats + i * num_floats; - - // copy row values - memcpy(dst, vals, stride); - - // calculate total cluster distance - sums[i] = sum_distance_squared(dst, set); - - // increment cluster count - i++; - } - } - - // FIXME: should probably be heap-allocated - int ints[num_clusters]; - memset(ints, 0, sizeof(ints)); - - // init cluster shape - const km_shape_t shape = { - .num_floats = num_floats, - .num_ints = 1, - }; - - // init cluster set - if (!km_set_init(cs, &shape, num_clusters)) { - // return failure - return false; - } - - // add data, return result - return km_set_push(cs, num_clusters, floats, ints); -} diff --git a/km-init-rand.c b/km-init-rand.c deleted file mode 100644 index 1f59da8..0000000 --- a/km-init-rand.c +++ /dev/null @@ -1,43 +0,0 @@ -#include // bool -#include // memset() -#include "util.h" -#include "km.h" - -// init a set with num_clusters clusters of shape num_floats by picking -// random cluster centers -bool -km_init_rand( - km_set_t * const cs, - const size_t num_clusters, - const km_set_t * const set, - km_rand_t * const rs -) { - // get number of floats from data set - const size_t num_floats = set->shape.num_floats; - - // init cluster shape - const km_shape_t shape = { - .num_floats = num_floats, - .num_ints = 1, - }; - - // generate random cluster centers - float floats[num_floats * num_clusters]; - if (!km_rand_get_floats(rs, num_floats * num_clusters, floats)) { - // return failure - return false; - } - - // FIXME: should probably be heap-allocated - int ints[num_clusters]; - memset(ints, 0, sizeof(ints)); - - // init cluster set - if (!km_set_init(cs, &shape, num_clusters)) { - // return failure - return false; - } - - // add data, return result - return km_set_push(cs, num_clusters, floats, ints); -} diff --git a/km-init.c b/km-init.c deleted file mode 100644 index 280ea81..0000000 --- a/km-init.c +++ /dev/null @@ -1,85 +0,0 @@ -#include // bool -#include // size_t -#include // strcmp() -#include "util.h" -#include "km.h" - -typedef bool (*km_init_fn_t)( - km_set_t *, - const size_t, - const km_set_t *, - km_rand_t * -); - -bool km_init_rand(km_set_t *, const size_t, const km_set_t *, km_rand_t *); -bool km_init_forgy(km_set_t *, const size_t, const km_set_t *, km_rand_t *); -bool km_init_kmeans(km_set_t *, const size_t, const km_set_t *, km_rand_t *); - -static const struct { - const km_init_type_t type; - const char * const name; - const km_init_fn_t init; -} TYPES[] = {{ - .name = "rand", - .type = KM_INIT_TYPE_RAND, - .init = km_init_rand, -}, { - .name = "forgy", - .type = KM_INIT_TYPE_FORGY, - .init = km_init_forgy, -}, { - .name = "kmeans", - .type = KM_INIT_TYPE_KMEANS, - .init = km_init_kmeans, -}}; - -#define NUM_TYPES (sizeof(TYPES) / sizeof(TYPES[0])) - -km_init_type_t -km_init_get_type( - const char * const s -) { - // find init method - for (size_t i = 0; i < NUM_TYPES; i++) { - if (!strcmp(s, TYPES[i].name)) { - // return type - return TYPES[i].type; - } - } - - // return failure - return KM_INIT_TYPE_LAST; -} - -static km_init_fn_t -get_init_func( - const km_init_type_t type -) { - // find init method - for (size_t i = 0; i < NUM_TYPES; i++) { - if (type == TYPES[i].type) { - // return init func - return TYPES[i].init; - } - } - - // return failure - return NULL; -} - -bool -km_init( - km_set_t * const cs, - const km_init_type_t init_type, - const size_t num_clusters, - const km_set_t * const set, - km_rand_t * const rs -) { - km_init_fn_t init_fn = get_init_func(init_type); - if (!init_fn) { - die("unknown cluster init method"); - } - - // call init, return result - return init_fn(cs, num_clusters, set, rs); -} diff --git a/km-load.c b/km-load.c deleted file mode 100644 index 3615aa9..0000000 --- a/km-load.c +++ /dev/null @@ -1,112 +0,0 @@ -#include // bool -#include // fscanf() -#include // strerror() -#include // errno -#include "util.h" -#include "km.h" - -#define FAIL(...) do { \ - if (cbs && cbs->on_error) { \ - char buf[1024]; \ - snprintf(buf, sizeof(buf), __VA_ARGS__); \ - cbs->on_error(buf, cb_data); \ - } \ - return false; \ -} while (0) - -_Bool -km_load( - FILE * const fh, - const km_load_cbs_t * const cbs, - void * const cb_data -) { - // read shape - km_shape_t shape = { 0, 0 }; - if (fscanf(fh, "%zu %zu", &(shape.num_floats), &(shape.num_ints)) != 2) { - FAIL("shape fscanf() failed: %s", strerror(errno)); - } - - if (cbs && cbs->on_shape) { - // emit shape - if (!cbs->on_shape(&shape, cb_data)) { - // return failure - return false; - } - } - - // alloc floats buffer - float *floats = NULL; - if (shape.num_floats > 0) { - floats = malloc(sizeof(float) * shape.num_floats); - if (!floats) { - FAIL("floats malloc() failed: %s", strerror(errno)); - } - } - - // alloc ints buffer - int *ints = NULL; - if (shape.num_ints > 0) { - ints = malloc(sizeof(int) * shape.num_ints); - if (!ints) { - FAIL("ints malloc() failed: %s", strerror(errno)); - } - } - - for (size_t row = 0; !feof(fh); row++) { - for (size_t i = 0; i < shape.num_floats; i++) { - if (fscanf(fh, " %f ", floats + i) != 1) { - FAIL("[%zu, %zu] float fscanf() failed: %s", row, i, strerror(errno)); - } - } - - // read ints - for (size_t i = 0; i < shape.num_ints; i++) { - if (fscanf(fh, " %d ", ints + i) != 1) { - FAIL("[%zu, %zu] int fscanf() failed: %s", row, i, strerror(errno)); - } - } - - if (cbs && cbs->on_row) { - // emit row - if (!cbs->on_row(floats, ints, cb_data)) { - // return failure - return false; - } - } - } - - if (shape.num_floats > 0) { - // free float buffer - free(floats); - } - - if (shape.num_ints > 0) { - // free int buffer - free(ints); - } - - // return success - return true; -} - -_Bool -km_load_path( - const char * const path, - const km_load_cbs_t * const cbs, - void * const cb_data -) { - // open file - FILE *fh = fopen(path, "rb"); - if (!fh) { - FAIL("fopen(\"%s\") failed: %s", path, strerror(errno)); - } - - // load file, get result - const bool r = km_load(fh, cbs, cb_data); - - // close file - fclose(fh); - - // return result - return r; -} diff --git a/km-rand-erand48.c b/km-rand-erand48.c deleted file mode 100644 index de84b93..0000000 --- a/km-rand-erand48.c +++ /dev/null @@ -1,89 +0,0 @@ -#include // bool -#define _DEFAULT_SOURCE -#include // drand48_r() -#include // fopen() -#include "util.h" -#include "km.h" - -// get N get_floats -static bool -on_get_floats( - km_rand_t * const rs, - const size_t num_vals, - float * const vals -) { - unsigned short * const data = rs->data; - - // generate results - for (size_t i = 0; i < num_vals; i++) { - vals[i] = erand48(data); - } - - // return success - return true; -} - -// fill sizes callback for system random source -static bool -on_get_sizes( - km_rand_t * const rs, - const size_t num_vals, - size_t * const vals -) { - unsigned short * const data = rs->data; - - // generate results - for (size_t i = 0; i < num_vals; i++) { - vals[i] = nrand48(data); - } - - // return success - return true; -} - -static void -on_fini( - km_rand_t * const rs -) { - unsigned short * data = rs->data; - - // free rng data - free(data); - - // clear pointer - rs->data = NULL; -} - -// path random source callbacks -static const km_rand_cbs_t -ERAND48_RAND_CBS = { - .get_floats = on_get_floats, - .get_sizes = on_get_sizes, - .fini = on_fini, -}; - -// init system random source (uses system rand()) -bool -km_rand_init_erand48( - km_rand_t * const rs, - const uint64_t seed -) { - // alloc rng data - unsigned short * const data = malloc(3 * sizeof(unsigned short)); - if (!data) { - // return failure - return false; - } - - // populate seed - data[0] = seed & 0xffff; - data[1] = (seed >> 16) & 0xffff; - data[2] = (seed >> 32) & 0xffff; - - // populate context - rs->cbs = &ERAND48_RAND_CBS; - rs->data = data; - - // return success - return true; -} diff --git a/km-rand-libc.c b/km-rand-libc.c deleted file mode 100644 index da1e4b0..0000000 --- a/km-rand-libc.c +++ /dev/null @@ -1,56 +0,0 @@ -#include // bool -#include "util.h" -#include "km.h" - -// libc random source get_floats -static bool -on_get_floats( - km_rand_t * const rs, - const size_t num_vals, - float * const vals -) { - UNUSED(rs); - - // generate results - for (size_t i = 0; i < num_vals; i++) { - vals[i] = 1.0 * rand() / RAND_MAX; - } - - // return success - return true; -} - -// fill sizes callback for system random source -static bool -on_get_sizes( - km_rand_t * const rs, - const size_t num_vals, - size_t * const vals -) { - UNUSED(rs); - - // generate results - for (size_t i = 0; i < num_vals; i++) { - vals[i] = rand(); - } - - // return success - return true; -} - -// system random source callbacks -static const km_rand_cbs_t -LIBC_RAND_CBS = { - .get_floats = on_get_floats, - .get_sizes = on_get_sizes, - .fini = NULL, -}; - -// init system random source (uses system rand()) -void -km_rand_init_libc( - km_rand_t * const rs -) { - rs->cbs = &LIBC_RAND_CBS; - rs->data = NULL; -} diff --git a/km-rand-path.c b/km-rand-path.c deleted file mode 100644 index 6597038..0000000 --- a/km-rand-path.c +++ /dev/null @@ -1,152 +0,0 @@ -#include // bool -#include // memset() -#include // fopen() -#include "util.h" -#include "km.h" - -#define MAX_BUF_ITEMS (1 << 20) - -typedef struct { - FILE *fh; - size_t ofs; - uint64_t buf[MAX_BUF_ITEMS]; -} ctx_t; - -static inline bool -refill( - ctx_t * const ctx -) { - ctx->ofs = 0; - // D("reading %zu bytes", sizeof(ctx->buf)); - return !!fread(ctx->buf, sizeof(ctx->buf), 1, ctx->fh); -} - -static bool -next( - ctx_t * const ctx, - uint64_t * const r -) { - if ((ctx->ofs >= MAX_BUF_ITEMS - 1) && !refill(ctx)) { - // return failure - return false; - } - - // get next value, increment offset - *r = ctx->buf[ctx->ofs]; - ctx->ofs++; - - // return success - return true; -} - -// get N get_floats -static bool -on_get_floats( - km_rand_t * const rs, - const size_t num_vals, - float * const vals -) { - ctx_t * const ctx = rs->data; - - // generate results - for (size_t i = 0; i < num_vals; i++) { - uint64_t v; - - if (!next(ctx, &v)) { - // return failure - return false; - } - - // add value - vals[i] = 1.0 * v / RAND_MAX; - } - - // return success - return true; -} - -// fill sizes callback for system random source -static bool -on_get_sizes( - km_rand_t * const rs, - const size_t num_vals, - size_t * const vals -) { - ctx_t * const ctx = rs->data; - - // generate results - for (size_t i = 0; i < num_vals; i++) { - uint64_t v; - - if (!next(ctx, &v)) { - // return failure - return false; - } - - // add value - vals[i] = v; - } - - // return success - return true; -} - -static void -on_fini( - km_rand_t * const rs -) { - ctx_t * const ctx = rs->data; - - // close file - fclose(ctx->fh); - - // free context - free(ctx); - - rs->data = NULL; -} - -// path random source callbacks -static const km_rand_cbs_t -PATH_RAND_CBS = { - .get_floats = on_get_floats, - .get_sizes = on_get_sizes, - .fini = on_fini, -}; - -// init system random source (uses system rand()) -bool -km_rand_init_path( - km_rand_t * const rs, - const char * const path -) { - // open file - FILE *fh = fopen(path, "rb"); - if (!fh) { - // return failure - return false; - } - - // alloc context - ctx_t * const ctx = malloc(sizeof(ctx_t)); - if (!ctx) { - // return failure - return false; - } - - // populate context - ctx->fh = fh; - - // fill buffer - if (!refill(ctx)) { - // return failure - return false; - } - - // populate random source - rs->cbs = &PATH_RAND_CBS; - rs->data = ctx; - - // return success - return true; -} diff --git a/km-rand.c b/km-rand.c deleted file mode 100644 index e777069..0000000 --- a/km-rand.c +++ /dev/null @@ -1,33 +0,0 @@ -#include -#include "util.h" -#include "km.h" - -// fill buffer with N random floats -bool -km_rand_get_floats( - km_rand_t * const rs, - const size_t num_floats, - float * const floats -) { - return rs->cbs->get_floats(rs, num_floats, floats); -} - -// fill buffer with N random size_ts -bool -km_rand_get_sizes( - km_rand_t * const rs, - const size_t num_sizes, - size_t * const sizes -) { - return rs->cbs->get_sizes(rs, num_sizes, sizes); -} - -// finalize random source -void -km_rand_fini( - km_rand_t * const rs -) { - if (rs->cbs->fini) { - rs->cbs->fini(rs); - } -} diff --git a/km-set.c b/km-set.c deleted file mode 100644 index 69ae6c0..0000000 --- a/km-set.c +++ /dev/null @@ -1,267 +0,0 @@ -#include // bool -#include // size_t -#include // memcpy() -#include // rand() -#include // errno -#include "util.h" -#include "km.h" - -#define MIN_ROWS (4096 / sizeof(float)) - -// grow data set -static bool -km_set_grow( - km_set_t * const set, - const size_t capacity -) { - float *floats = NULL; - const size_t floats_size = sizeof(float) * set->shape.num_floats * capacity; - // fprintf(stderr, "floats_size = %zu\n", floats_size); - if (floats_size > 0) { - // alloc floats - floats = realloc(set->floats, floats_size); - if (!floats) { - // return failure - return false; - } - } - - int *ints = NULL; - const size_t ints_size = sizeof(int) * set->shape.num_ints * capacity; - // fprintf(stderr, "ints_size = %zu\n", ints_size); - if (ints_size > 0) { - // alloc ints - ints = realloc(set->ints, ints_size); - if (!ints) { - // return failure - return false; - } - } - - // update set - set->floats = floats; - set->ints = ints; - set->capacity = capacity; - - // return success - return true; -} - -// init data set with shape and initial size -bool -km_set_init( - km_set_t * const set, - const km_shape_t * const shape, - const size_t row_capacity -) { - // alloc bounds - float * const bounds = malloc(2 * sizeof(float) * shape->num_floats); - if (!bounds) { - // return error - return false; - } - - set->state = KM_SET_STATE_INIT; - set->floats = NULL; - set->ints = NULL; - set->shape = *shape; - set->num_rows = 0; - set->capacity = 0; - set->bounds = bounds; - - return km_set_grow(set, MAX(MIN_ROWS, row_capacity + 1)); -} - -// finalize data set -void -km_set_fini(km_set_t * const set) { - if (set->state == KM_SET_STATE_FINI) { - return; - } - - if (set->bounds) { - // free bounds - free(set->bounds); - set->bounds = NULL; - } - - if (set->floats) { - // free floats - free(set->floats); - set->floats = NULL; - } - - if (set->ints) { - // free ints - free(set->ints); - set->ints = NULL; - } - - // shrink capacity - set->capacity = 0; - - // set state - set->state = KM_SET_STATE_FINI; -} - -// append rows to data set, growing set if necessary -bool -km_set_push( - km_set_t * const set, - const size_t num_rows, - const float * const floats, - const int * const ints -) { - // check state - if (set->state != KM_SET_STATE_INIT) { - // return failure - return false; - } - - const size_t capacity_needed = set->num_rows + num_rows; - // FIXME: potential overflow here - if (capacity_needed >= set->capacity) { - // crappy growth algorithm - const size_t new_capacity = 2 * capacity_needed + 1; - - // resize storage - if (!km_set_grow(set, MAX(MIN_ROWS, new_capacity))) { - return false; - } - } - - // copy floats - const size_t num_floats = set->shape.num_floats; - if (num_floats > 0) { - float * const dst = set->floats + num_floats * set->num_rows; - const size_t stride = sizeof(float) * num_floats; - - // copy floats - memcpy(dst, floats, stride * num_rows); - - if (!set->num_rows) { - // there were no rows, so populate bounds with first row - memcpy(set->bounds, floats, stride); - memcpy(set->bounds + num_floats, floats, stride); - } - - for (size_t i = 0; i < num_rows; i++) { - for (size_t j = 0; j < num_floats; j++) { - const float val = floats[i * num_floats + j]; - - if (val < set->bounds[j]) { - // update min bound - set->bounds[j] = val; - } - - if (val > set->bounds[num_floats + j]) { - // update max bound - set->bounds[num_floats + j] = val; - } - } - } - } - - // copy ints - const size_t num_ints = set->shape.num_ints; - if (num_ints > 0) { - int * const dst = set->ints + num_ints * set->num_rows; - const size_t stride = sizeof(int) * num_ints; - - // copy ints - memcpy(dst, ints, stride * num_rows); - } - - // increment row count - set->num_rows += num_rows; - - // return success - return true; -} - -bool -km_set_copy( - km_set_t * const dst, - const km_set_t * const src -) { - if (src->state != KM_SET_STATE_INIT && src->state != KM_SET_STATE_NORMALIZED) { - // return failure - D("invalid state"); - return false; - } - - // init dst set - if (!km_set_init(dst, &(src->shape), src->num_rows)) { - // return failure - D("km_set_init()"); - return false; - } - - // copy floats - if (src->shape.num_floats > 0) { - const size_t stride = sizeof(float) * src->shape.num_floats; - - // copy floats - memcpy(dst->floats, src->floats, stride * src->num_rows); - - // copy bounds - memcpy(dst->bounds, src->bounds, 2 * stride); - } - - // copy ints - const size_t num_ints = src->shape.num_ints; - if (num_ints > 0) { - const size_t stride = sizeof(int) * num_ints; - - // copy ints - memcpy(dst->ints, src->ints, stride * src->num_rows); - } - - // increment row count - dst->num_rows = src->num_rows; - - // return success - return true; -} - -bool -km_set_normalize( - km_set_t * const set -) { - const size_t num_floats = set->shape.num_floats; - - // check set state - if (set->state != KM_SET_STATE_INIT) { - // return failure - return false; - } - - // normalize values - for (size_t i = 0; i < set->num_rows; i++) { - for (size_t j = 0; j < num_floats; j++) { - const size_t ofs = i * num_floats + j; - const float val = set->floats[ofs], - min = set->bounds[j], - max = set->bounds[num_floats + j]; - - // normalize and write value - set->floats[ofs] = (val - min) / (max - min); - } - } - - // set state - set->state = KM_SET_STATE_NORMALIZED; - - // return success - return true; -} - -// get row from data set -float * -km_set_get_row( - const km_set_t * const set, - const size_t i -) { - const size_t num_floats = set->shape.num_floats; - return set->floats + i * num_floats; -} diff --git a/km-solve.c b/km-solve.c deleted file mode 100644 index f579a41..0000000 --- a/km-solve.c +++ /dev/null @@ -1,200 +0,0 @@ -#include // bool -#include // memset() -#include // FLT_MAX -#include // sqrt() -#include "util.h" -#include "km.h" - -#define MIN_CLUSTER_DISTANCE 0.00001 - -// alloc and initialize row map -static km_row_map_t * -km_row_map_init( - const size_t num_rows -) { - // alloc row map - km_row_map_t * const row_map = malloc(sizeof(km_row_map_t) * num_rows); - - // check for error - if (!row_map) { - // return failure - return false; - } - - // init row map - for (size_t i = 0; i < num_rows; i++) { - // setting distances to maximum - row_map[i].d2 = FLT_MAX; - row_map[i].d2_near = FLT_MAX; - } - - // return row map - return row_map; -} - -static void -km_row_map_fini( - km_row_map_t * const row_map -) { - free(row_map); -} - -// use k-means to iteratively update the cluster centroids until there -// are no more changes to the centroids -bool -km_solve( - km_set_t * const cs, - const km_set_t * const set, - const km_solve_cbs_t * const cbs, - void * const cb_data -) { - const size_t num_clusters = cs->num_rows, - num_floats = set->shape.num_floats; - - // row map: row => distance, cluster ID - km_row_map_t *row_map = km_row_map_init(set->num_rows); - if (!row_map) { - // return failure - return false; - } - - // calculate clusters by doing the following: - // * walk all clusters and all rows - // * if we find a closer cluster, move row to cluster - // * if there were changes to any cluster, then calculate a new - // centroid for each cluster by averaging the cluster rows - // * repeat until there are no more changes - bool changed = false; - do { - // no changes yet - changed = false; - - for (size_t i = 0; i < num_clusters; i++) { - // get the floats for this cluster - const float * const floats = km_set_get_row(cs, i); - - for (size_t j = 0; j < set->num_rows; j++) { - // get row values - const float * const vals = km_set_get_row(set, j); - - // calculate the distance squared between row and cluster - const float d2 = distance_squared(num_floats, floats, vals); - - if (d2 < row_map[j].d2) { - // row is closer to this cluster, update row map - row_map[j].d2 = d2; - row_map[j].cluster = i; - - // flag change - changed = true; - } - - if ((row_map[j].cluster != i) && (d2 < row_map[j].d2_near)) { - row_map[j].d2_near = d2; - - // flag change - changed = true; - } - } - } - - if (changed) { - // if there were changes, then we need to calculate the new - // cluster centers - - // calculate new center - for (size_t i = 0; i < num_clusters; i++) { - size_t num_rows = 0; - float * const floats = km_set_get_row(cs, i); - memset(floats, 0, sizeof(float) * num_floats); - - for (size_t j = 0; j < set->num_rows; j++) { - const float * const row_floats = km_set_get_row(set, j); - - if (row_map[j].cluster == i) { - // calculate numerator for average - for (size_t k = 0; k < num_floats; k++) { - floats[k] += row_floats[k]; - } - - // increment denominator for average - num_rows++; - } - } - - // save number of rows in this cluster - cs->ints[i] = num_rows; - - if (num_rows > 0) { - for (size_t k = 0; k < num_floats; k++) { - // divide by denominator to get average - floats[k] /= num_rows; - } - } - } - } - - if (cbs && cbs->on_step) { - // pass clusters and row map to step callback - cbs->on_step(cs, row_map, cb_data); - } - } while (changed); - - if (cbs && cbs->on_stats) { - float sum = 0, - silouette = 0, - mean_dists[num_clusters], - mean_nears[num_clusters]; - - memset(mean_dists, 0, sizeof(mean_dists)); - memset(mean_nears, 0, sizeof(mean_nears)); - - // calculate sum of distances across all clusters - for (size_t i = 0; i < set->num_rows; i++) { - sum += row_map[i].d2; - } - - // calculate mean numerators and silouette - for (size_t i = 0; i < set->num_rows; i++) { - // distance squared (d2) to center of this cluster - mean_dists[row_map[i].cluster] += row_map[i].d2; - - // distance squared (d2) to center of nearest cluster - mean_nears[row_map[i].cluster] += row_map[i].d2_near; - - // calculate silouette denominator - // (https://en.wikipedia.org/wiki/Silhouette_%28clustering%29) - const float delta = row_map[i].d2_near - row_map[i].d2; - if (fabsf(delta) > MIN_CLUSTER_DISTANCE) { - silouette += delta / MAX(row_map[i].d2, row_map[i].d2_near); - } - } - - // finalize means (divide by row count) - for (size_t i = 0; i < num_clusters; i++) { - mean_dists[i] = (cs->ints[i]) ? (sqrt(mean_dists[i]) / cs->ints[i]) : 0; - mean_nears[i] = (cs->ints[i]) ? (sqrt(mean_nears[i]) / cs->ints[i]) : 0; - } - - // finalize silouette - silouette /= set->num_rows; - - // build stats - const km_solve_stats_t stats = { - .sum = sum, - .silouette = silouette, - .mean_dists = mean_dists, - .mean_nears = mean_nears, - .num_clusters = num_clusters, - }; - - // emit means - cbs->on_stats(cs, &stats, cb_data); - } - - // free row_map - km_row_map_fini(row_map); - - // return success - return true; -} diff --git a/km.h b/km.h deleted file mode 100644 index 5a34261..0000000 --- a/km.h +++ /dev/null @@ -1,230 +0,0 @@ -#ifndef KM_H -#define KM_H - -#include // size_t -#include // uint8_t -#include // FILE - -// forward typedef for callbacks -typedef struct km_rand_t_ km_rand_t; - -// random number source callbacks -typedef struct { - _Bool (*get_floats)(km_rand_t * const, const size_t, float * const); - _Bool (*get_sizes)(km_rand_t * const, const size_t, size_t * const); - void (*fini)(km_rand_t * const); -} km_rand_cbs_t; - -// random number source -struct km_rand_t_ { - const km_rand_cbs_t *cbs; - void *data; -}; - -// fill buffer with N random floats -_Bool km_rand_get_floats(km_rand_t * const, const size_t, float * const); - -// fill buffer with N random size_ts -_Bool km_rand_get_sizes(km_rand_t * const, const size_t, size_t * const); - -// finalize random source -void km_rand_fini(km_rand_t * const); - -// init libc rng (use libc rand()) -void km_rand_init_libc(km_rand_t *); - -// init file source rng (e.g. "/dev/urandom") -_Bool km_rand_init_path(km_rand_t * const, const char *); - -// init erand48 rng (use erand48()) -_Bool km_rand_init_erand48(km_rand_t * const, const uint64_t); - -// shape of data set -typedef struct { - size_t num_floats, - num_ints; -} km_shape_t; - -typedef enum { - KM_SET_STATE_INIT, - KM_SET_STATE_NORMALIZED, - KM_SET_STATE_FINI, - KM_SET_STATE_LAST, -} km_set_state_t; - -// data set -typedef struct { - km_set_state_t state; - - float *floats, - *bounds; - int *ints; - - km_shape_t shape; - - size_t num_rows, - capacity; -} km_set_t; - -// init data set with shape and initial size -_Bool km_set_init(km_set_t * const, const km_shape_t *, const size_t); -// finalize data set -void km_set_fini(km_set_t * const); - -// append rows to data set, growing set if necessary -_Bool km_set_push(km_set_t *, const size_t, const float *, const int *); - -// deep copy data set -_Bool km_set_copy(km_set_t * const, const km_set_t * const); - -// normalize data set -_Bool km_set_normalize(km_set_t * const); - -// get pointer to data set row -float *km_set_get_row(const km_set_t *, const size_t); - -typedef struct { - float d2, - d2_near; - size_t cluster; -} km_row_map_t; - -typedef struct { - const float sum; - - const float silouette; - - const float *mean_dists; - const float *mean_nears; - - const size_t num_clusters; -} km_solve_stats_t; - -typedef struct { - void (*on_step)(const km_set_t *, const km_row_map_t *, void *); - void (*on_stats)(const km_set_t *, const km_solve_stats_t *, void *); -} km_solve_cbs_t; - -// use k-means to iteratively update the cluster centroids until there -// are no more changes to the centroids -_Bool -km_solve( - km_set_t * const, - const km_set_t * const, - const km_solve_cbs_t * const, - void * const -); - -typedef struct { - const km_set_t * const cluster_set; - - float distance_sum, - silouette, - mean_cluster_size; - - const size_t num_clusters, - num_empty_clusters; -} km_find_data_t; - -typedef _Bool (*km_find_init_cb_t)( - km_set_t * const cluster_set, - const size_t num_clusters, - const km_set_t * const data_set, - void *cb_data -); - -typedef _Bool (*km_find_fini_cb_t)( - km_set_t * const cluster_set, - void *cb_data -); - -typedef _Bool (*km_find_best_cb_t)( - const float score, - const km_set_t * const cluster_set, - void *cb_data -); - -typedef void (*km_find_data_cb_t)( - const km_find_data_t * const, - void * -); - -typedef struct { - // maximum number of clusters - size_t max_clusters; - - // maximum number of tests per cluster - size_t num_tests; - - km_find_init_cb_t on_init; - km_find_fini_cb_t on_fini; - km_find_data_cb_t on_data; - km_find_best_cb_t on_best; -} km_find_cbs_t; - -// repeatedly test different cluster sizes and report results -_Bool -km_find( - const km_set_t * const set, - const km_find_cbs_t * const cbs, - void * const cb_data -); - -// render set to rgb buffer -void -km_set_draw( - const km_set_t * const set, - uint8_t * const rgb, - const size_t width, - const size_t height, - const int dot_size, - const uint32_t color -); - -// load callbacks -typedef struct { - _Bool (*on_shape)(const km_shape_t * const, void *); - _Bool (*on_row)(const float *, const int *, void *); - void (*on_error)(const char * const, void *); -} km_load_cbs_t; - -// load set from file handle -_Bool -km_load( - FILE * const fh, - const km_load_cbs_t * const cbs, - void * const cb_data -); - -// load set from path -_Bool km_load_path( - const char * const path, - const km_load_cbs_t * const cbs, - void * const cb_data -); - -typedef enum { - // rand: init cluster set by picking random centroids - KM_INIT_TYPE_RAND, - - // forgy: init cluster set by picking random rows from set - KM_INIT_TYPE_FORGY, - - // kmeans++: ... - KM_INIT_TYPE_KMEANS, - - KM_INIT_TYPE_LAST, -} km_init_type_t; - -km_init_type_t km_init_get_type(const char * const); - -// init cluster set using given method -_Bool km_init( - km_set_t * const, - const km_init_type_t, - const size_t, - const km_set_t * const, - km_rand_t * const -); - -#endif /* KM_H */ diff --git a/main.c b/main.c deleted file mode 100644 index 55056b7..0000000 --- a/main.c +++ /dev/null @@ -1,398 +0,0 @@ -#include // bool -#include // fprintf() -#include // EXIT_{FAILURE,SUCCESS} -#include // exit() -#include // memset() -#include // fabsf() - -#define STB_IMAGE_WRITE_IMPLEMENTATION -#define STB_ONLY_PNG -#include "stb_image_write.h" - -#include "util.h" -#include "km.h" - -#define MAX_CLUSTERS 10 -#define NUM_TESTS 300 -#define MAX_BEST 10 - -#define IM_WIDTH 128 -#define IM_HEIGHT 128 -#define IM_STRIDE (3 * IM_WIDTH) - -typedef struct { - float score; - km_set_t set; -} best_item_t; - -typedef struct { - // cluster initialization method - km_init_type_t init_type; - - // random number source - km_rand_t rs; - - struct { - float distance, - silouette, - cluster_size; - size_t num_empty; - } rows[MAX_CLUSTERS - 2]; - - // best clusters - best_item_t best[MAX_BEST]; - size_t num_best; -} ctx_t; - -static int -best_score_cmp( - const void * const ap, - const void * const bp -) { - const best_item_t * const a = ap; - const best_item_t * const b = bp; - - return (a->score > b->score) ? 1 : -1; -} - -static void -ctx_best_sort( - ctx_t * const ctx -) { - // sort best sets by ascending score (worst to best) - qsort( - ctx->best, - MIN(ctx->num_best, MAX_BEST), - sizeof(best_item_t), - best_score_cmp - ); -} - -static void -ctx_best_walk( - const ctx_t * const ctx, - void (*on_best)(const km_set_t * const, const size_t, const float, void *), - void * const cb_data -) { - if (on_best) { - // walk best sets and emit each one - for (size_t i = 0; i < MIN(ctx->num_best, MAX_BEST); i++) { - on_best(&(ctx->best[i].set), i, ctx->best[i].score, cb_data); - } - } -} - -static bool -load_on_shape( - const km_shape_t * const shape, - void * const cb_data -) { - km_set_t * const set = cb_data; - - // D("shape: %zu floats, %zu ints", shape->num_floats, shape->num_ints); - - // init set - if (!km_set_init(set, shape, 100)) { - die("km_set_init() failed"); - } - - // return success - return true; -} - -static bool -load_on_row( - const float * const floats, - const int * const ints, - void * const cb_data -) { - km_set_t * const set = cb_data; - - // push row - if (!km_set_push(set, 1, floats, ints)) { - die("km_set_push_rows() failed"); - } - - // return success - return true; -} - -static void -load_on_error( - const char * const err, - void * const cb_data -) { - UNUSED(cb_data); - die("load failed: %s", err); -} - -static const km_load_cbs_t -LOAD_CBS = { - .on_shape = load_on_shape, - .on_row = load_on_row, - .on_error = load_on_error, -}; - -static bool -find_on_init( - km_set_t * const cs, - const size_t num_clusters, - const km_set_t * const set, - void *cb_data -) { - ctx_t * const ctx = cb_data; - return km_init(cs, ctx->init_type, num_clusters, set, &(ctx->rs)); -} - -static bool -find_on_fini( - km_set_t * const cs, - void *cb_data -) { - UNUSED(cb_data); - km_set_fini(cs); - return true; -} - -static void -find_on_data( - const km_find_data_t * const data, - void *cb_data -) { - ctx_t * const ctx = cb_data; - const size_t ofs = data->num_clusters - 2; - - ctx->rows[ofs].distance += data->distance_sum; - ctx->rows[ofs].silouette += data->silouette; - ctx->rows[ofs].cluster_size += data->mean_cluster_size; - ctx->rows[ofs].num_empty += data->num_empty_clusters; -} - -static bool -find_on_best( - const float score, - const km_set_t * const cs, - void *cb_data -) { - ctx_t * const ctx = cb_data; - - D("new best: score = %0.3f, num_clusters = %zu", score, cs->num_rows); - - // get pointer to destination set - // (note: data->best is a ring buffer) - const size_t ofs = ctx->num_best % MAX_BEST; - ctx->best[ofs].score = score; - km_set_t * const dst = &(ctx->best[ofs].set); - - if (ctx->num_best >= MAX_BEST) { - // finalize old best data set - km_set_fini(dst); - } - - // copy data set to best ring buffer - if (!km_set_copy(dst, cs)) { - die("km_set_copy()"); - } - - // increment best count - ctx->num_best++; - - // return success - return true; -} - -// init find config -static const km_find_cbs_t -FIND_CBS = { - .max_clusters = MAX_CLUSTERS, - .num_tests = NUM_TESTS, - - .on_init = find_on_init, - .on_fini = find_on_fini, - .on_data = find_on_data, - .on_best = find_on_best, -}; - -static void -ctx_csv_print_row( - const ctx_t * const ctx, - FILE * const fh, - const size_t i -) { - const size_t num_clusters = i + 2; - const float mean_distance = ctx->rows[i].distance / NUM_TESTS, - mean_cluster_size = ctx->rows[i].cluster_size / NUM_TESTS, - mean_empty = 1.0 * ctx->rows[i].num_empty / NUM_TESTS, - score = ctx->rows[i].silouette / NUM_TESTS; - - // print result - fprintf(fh, "%zu,%0.3f,%0.3f,%0.3f,%0.3f\n", - num_clusters, - score, - mean_distance, - mean_cluster_size, - mean_empty - ); -} - -static void -ctx_csv_print( - const ctx_t * const ctx, - FILE * const fh -) { - // print headers - fprintf(fh, "#,score,distance,size,empty\n"); - - // print rows - for (size_t i = 0; i < MAX_CLUSTERS - 2; i++) { - ctx_csv_print_row(ctx, fh, i); - } -} - -// static image data buffer -static uint8_t im_data[3 * IM_WIDTH * IM_HEIGHT]; - -static void -save_on_best( - const km_set_t * const set, - const size_t rank, - const float score, - void * const cb_data -) { - const ctx_t * const ctx = cb_data; - UNUSED(score); - - // convert rank to channel brightness - const uint8_t ch = 0x33 + (0xff - 0x33) * (1.0 * rank + 1) / (MAX_BEST); - const uint8_t shift = (rank == MIN(ctx->num_best, MAX_BEST) - 1) ? 8 : 16; - const uint32_t color = (ch & 0xff) << shift; - // const uint32_t color = 0xff0000; - // D("rank = %zu, score = %0.3f, size = %zu, color = %06x", rank, score, set->num_rows, color); - - // draw clusters - km_set_draw(set, im_data, IM_WIDTH, IM_HEIGHT, 3, color); -} - -static void -ctx_save_png( - const ctx_t * const ctx, - const char * const png_path, - const km_set_t * const set -) { - // clear image data to white - memset(im_data, 0xff, sizeof(im_data)); - - // draw data points - km_set_draw(set, im_data, IM_WIDTH, IM_HEIGHT, 1, 0x000000); - if (!stbi_write_png(png_path, IM_WIDTH, IM_HEIGHT, 3, im_data, IM_STRIDE)) { - die("stbi_write_png(\"%s\")", png_path); - } - - // draw best cluster points - ctx_best_walk(ctx, save_on_best, (void*) ctx); - - // save png - if (!stbi_write_png(png_path, IM_WIDTH, IM_HEIGHT, 3, im_data, IM_STRIDE)) { - die("stbi_write_png(\"%s\")", png_path); - } -} - -static void -ctx_best_print_on_best( - const km_set_t * const set, - const size_t rank, - const float score, - void * const cb_data -) { - FILE * const fh = cb_data; - - fprintf(fh, "rank = %zu, score = %0.3f, num_clusters = %zu: [\n", rank, score, set->num_rows); - for (size_t i = 0; i < set->num_rows; i++) { - const float * const vals = km_set_get_row(set, i); - fprintf(fh, " ["); - - for (size_t j = 0; j < set->shape.num_floats; j++) { - fprintf(fh, "%s%0.3f", (j > 0) ? ", " : "", vals[j]); - } - - fprintf(fh, "], (%d rows)\n", set->ints[i]); - } - fprintf(fh, "]\n"); -} - -static void -ctx_best_print( - const ctx_t * const ctx, - FILE * const fh -) { - ctx_best_walk(ctx, ctx_best_print_on_best, fh); -} - -static const char USAGE_FORMAT[] = - "Usage: %s [init] [data_path] \n" - "\n" - "Arguments:\n" - "* init: Cluster init method (one of \"rand\" or \"set\").\n" - "* data_path: Path to input data file.\n" - "* png_path: Path to output file (optional).\n" - ""; - -int main(int argc, char *argv[]) { - // check command-line arguments - if (argc < 3) { - fprintf(stderr, USAGE_FORMAT, argv[0]); - return EXIT_FAILURE; - } - - // get command-line arguments - const char * const init_type_name = argv[1]; - const char * const data_path = argv[2]; - const char * const png_path = (argc > 3) ? argv[3] : NULL; - - // init random seed - srand(getpid()); - - // init context - ctx_t ctx; - memset(&ctx, 0, sizeof(ctx_t)); - ctx.init_type = km_init_get_type(init_type_name); - - // init ctx rng - km_rand_init_erand48(&(ctx.rs), rand()); - - // init data set - km_set_t set; - if (!km_load_path(data_path, &LOAD_CBS, &set)) { - die("km_load_path(\"%s\") failed", data_path); - } - - // normalize data set - if (!km_set_normalize(&set)) { - die("km_set_normalize() failed"); - } - - // find best solutions - if (!km_find(&set, &FIND_CBS, &ctx)) { - die("km_find()"); - } - - // print csv - ctx_csv_print(&ctx, stdout); - - // sort best results from lowest to highest - ctx_best_sort(&ctx); - - // print best - ctx_best_print(&ctx, stdout); - - if (png_path) { - // save png of normalized data set and best clusters - ctx_save_png(&ctx, png_path, &set); - } - - // finalize data set - km_set_fini(&set); - - // return success - return 0; -} diff --git a/src/km-draw.c b/src/km-draw.c new file mode 100644 index 0000000..ddf1b26 --- /dev/null +++ b/src/km-draw.c @@ -0,0 +1,31 @@ +#include // size_t +#include "util.h" +#include "km.h" + +void +km_set_draw( + const km_set_t * const set, + uint8_t * const rgb, + const size_t width, + const size_t height, + const int dot_size, + const uint32_t color +) { + for (size_t i = 0; i < set->num_rows; i++) { + const float *row = km_set_get_row(set, i); + + for (int yo = 0; yo < dot_size; yo++) { + for (int xo = 0; xo < dot_size; xo++) { + const int x = (width - 1) * row[0] - dot_size / 2 + xo, + y = (height - 1) * row[1] - dot_size / 2 + yo, + ofs = 3 * (width * y + x); + + if (x >= 0 && x < (int) width && y >= 0 && y < (int) height) { + rgb[ofs + 0] = (color & 0xff0000) >> 16; + rgb[ofs + 1] = (color & 0x00ff00) >> 8; + rgb[ofs + 2] = (color & 0x0000ff); + } + } + } + } +} diff --git a/src/km-find.c b/src/km-find.c new file mode 100644 index 0000000..657c645 --- /dev/null +++ b/src/km-find.c @@ -0,0 +1,115 @@ +#include // bool +#include // fabsf() +#include "util.h" +#include "km.h" + +#define MIN_CLUSTER_DISTANCE 0.0001 + +static float +get_mean_cluster_size( + const km_set_t * const set +) { + float sum = 0; + size_t num_filled = 0; + + for (size_t i = 0; i < set->num_rows; i++) { + if (set->ints[i] > 0) { + sum += set->ints[i]; + num_filled++; + } + } + + return (num_filled > 0) ? (sum / num_filled) : 0; +} + +static void +find_solve_on_stats( + const km_set_t * const set, + const km_solve_stats_t * const stats, + void * const cb_data +) { + km_find_data_t * const find_data = cb_data; + UNUSED(set); + + // save total sum and silouette + find_data->distance_sum = stats->sum; + find_data->silouette = stats->silouette; +} + +static const km_solve_cbs_t +FIND_SOLVE_CBS = { + .on_stats = find_solve_on_stats, +}; + +bool +km_find( + const km_set_t * const set, + const km_find_cbs_t * const cbs, + void * const cb_data +) { + // check init callback + if (!cbs->on_init) { + return false; + } + + // check fini callback + if (!cbs->on_fini) { + return false; + } + + // check data callback + if (!cbs->on_data) { + return false; + } + + float best_silouette = -2.0; + for (size_t i = 2; i < cbs->max_clusters; i++) { + for (size_t j = 0; j < cbs->num_tests; j++) { + // init cluster set + km_set_t cs; + if (!cbs->on_init(&cs, i, set, cb_data)) { + // return failure + return false; + } + + // init find data + km_find_data_t find_data = { + .cluster_set = &cs, + .num_clusters = i, + }; + + // solve test + // (populates sum and silouette) + if (!km_solve(&cs, set, &FIND_SOLVE_CBS, &find_data)) { + // return failure + return false; + } + + // populate mean cluster size + find_data.mean_cluster_size = get_mean_cluster_size(&cs); + + // emit result + cbs->on_data(&find_data, cb_data); + + if (find_data.silouette > best_silouette) { + // emit new best result + if (cbs->on_best && !cbs->on_best(find_data.silouette, &cs, cb_data)) { + // return failure + return false; + } + + // update best silouette + best_silouette = find_data.silouette; + } + + // finalize cluster set + if (!cbs->on_fini(&cs, cb_data)) { + // return failure + return false; + } + } + } + + // return success + return true; +} diff --git a/src/km-init-forgy.c b/src/km-init-forgy.c new file mode 100644 index 0000000..5ff0c12 --- /dev/null +++ b/src/km-init-forgy.c @@ -0,0 +1,51 @@ +#include // bool +#include // memset() +#include "util.h" +#include "km.h" + +// init a set with num_clusters clusters of shape num_floats by picking +// random initial points from the set +bool +km_init_forgy( + km_set_t * const cs, + const size_t num_clusters, + const km_set_t * const set, + km_rand_t * const rs +) { + const size_t num_floats = set->shape.num_floats, + stride = sizeof(float) * num_floats; + + // init cluster shape + const km_shape_t shape = { + .num_floats = num_floats, + .num_ints = 1, + }; + + // get random row offsets + size_t rows[num_clusters]; + if (!km_rand_get_sizes(rs, num_clusters, rows)) { + // return failure + return false; + } + + // generate random cluster centers + float floats[num_floats * num_clusters]; + for (size_t i = 0; i < num_clusters; i++) { + const size_t row_num = rows[i] % set->num_rows; + const float * const row_floats = km_set_get_row(set, row_num); + memcpy(floats + i * num_floats, row_floats, stride); + } + + // FIXME: should probably be heap-allocated + int ints[num_clusters]; + memset(ints, 0, sizeof(ints)); + + // init cluster set + if (!km_set_init(cs, &shape, num_clusters)) { + // return failure + return false; + } + + // add data, return result + return km_set_push(cs, num_clusters, floats, ints); +} diff --git a/src/km-init-kmeans.c b/src/km-init-kmeans.c new file mode 100644 index 0000000..f0c1b5f --- /dev/null +++ b/src/km-init-kmeans.c @@ -0,0 +1,124 @@ +#include // bool +#include // memset() +#include // FLT_MAX +#include "util.h" +#include "km.h" + +// sum the squared distance of every row in set from this point +static inline float +sum_distance_squared( + const float * const floats, + const km_set_t * const set +) { + float r = 0; + + // sum squared distances + for (size_t i = 0; i < set->num_rows; i++) { + const float * const vals = km_set_get_row(set, i); + r += distance_squared(set->shape.num_floats, floats, vals); + } + + // return result + return r; +} + +static const float * +get_random_row( + const km_set_t * const set, + km_rand_t * const rs +) { + // get random offset + size_t ofs = 0; + if (!km_rand_get_sizes(rs, 1, &ofs)) { + die("km_rand_fill_sizes()"); + } + + return km_set_get_row(set, ofs % set->num_rows); +} + +// init a cluster set of num_clusters by picking the first cluster from +// the set of points at random +// random initial points from the set +bool +km_init_kmeans( + km_set_t * const cs, + const size_t num_clusters, + const km_set_t * const set, + km_rand_t * const rs +) { + const size_t num_floats = set->shape.num_floats, + stride = sizeof(float) * num_floats; + + // row values (filled below) + float sums[num_clusters]; + float floats[num_floats * num_clusters]; + + // pick first row randomly + { + // get random row + const float * const vals = get_random_row(set, rs); + + // copy row values to floats buffer + memcpy(floats, vals, stride); + + // copy row values to floats buffer + sums[0] = sum_distance_squared(floats, set); + } + + for (size_t i = 1; i < num_clusters;) { + // get random row + const float * const vals = get_random_row(set, rs); + + // calculate squared distance to nearest cluster + size_t ofs = 0; + float min_d2 = FLT_MAX; + for (size_t j = 0; j < i; j++) { + const float d2 = distance_squared(num_floats, vals, floats + j * stride); + if (d2 < min_d2) { + ofs = j; + min_d2 = d2; + } + } + + // get a random floating point value + float rand_val = 0; + if (!km_rand_get_floats(rs, 1, &rand_val)) { + // return failure + return false; + } + + // check random value + if (rand_val / sums[ofs] < min_d2) { + // get destination + float * const dst = floats + i * num_floats; + + // copy row values + memcpy(dst, vals, stride); + + // calculate total cluster distance + sums[i] = sum_distance_squared(dst, set); + + // increment cluster count + i++; + } + } + + // FIXME: should probably be heap-allocated + int ints[num_clusters]; + memset(ints, 0, sizeof(ints)); + + // init cluster shape + const km_shape_t shape = { + .num_floats = num_floats, + .num_ints = 1, + }; + + // init cluster set + if (!km_set_init(cs, &shape, num_clusters)) { + // return failure + return false; + } + + // add data, return result + return km_set_push(cs, num_clusters, floats, ints); +} diff --git a/src/km-init-rand.c b/src/km-init-rand.c new file mode 100644 index 0000000..1f59da8 --- /dev/null +++ b/src/km-init-rand.c @@ -0,0 +1,43 @@ +#include // bool +#include // memset() +#include "util.h" +#include "km.h" + +// init a set with num_clusters clusters of shape num_floats by picking +// random cluster centers +bool +km_init_rand( + km_set_t * const cs, + const size_t num_clusters, + const km_set_t * const set, + km_rand_t * const rs +) { + // get number of floats from data set + const size_t num_floats = set->shape.num_floats; + + // init cluster shape + const km_shape_t shape = { + .num_floats = num_floats, + .num_ints = 1, + }; + + // generate random cluster centers + float floats[num_floats * num_clusters]; + if (!km_rand_get_floats(rs, num_floats * num_clusters, floats)) { + // return failure + return false; + } + + // FIXME: should probably be heap-allocated + int ints[num_clusters]; + memset(ints, 0, sizeof(ints)); + + // init cluster set + if (!km_set_init(cs, &shape, num_clusters)) { + // return failure + return false; + } + + // add data, return result + return km_set_push(cs, num_clusters, floats, ints); +} diff --git a/src/km-init.c b/src/km-init.c new file mode 100644 index 0000000..280ea81 --- /dev/null +++ b/src/km-init.c @@ -0,0 +1,85 @@ +#include // bool +#include // size_t +#include // strcmp() +#include "util.h" +#include "km.h" + +typedef bool (*km_init_fn_t)( + km_set_t *, + const size_t, + const km_set_t *, + km_rand_t * +); + +bool km_init_rand(km_set_t *, const size_t, const km_set_t *, km_rand_t *); +bool km_init_forgy(km_set_t *, const size_t, const km_set_t *, km_rand_t *); +bool km_init_kmeans(km_set_t *, const size_t, const km_set_t *, km_rand_t *); + +static const struct { + const km_init_type_t type; + const char * const name; + const km_init_fn_t init; +} TYPES[] = {{ + .name = "rand", + .type = KM_INIT_TYPE_RAND, + .init = km_init_rand, +}, { + .name = "forgy", + .type = KM_INIT_TYPE_FORGY, + .init = km_init_forgy, +}, { + .name = "kmeans", + .type = KM_INIT_TYPE_KMEANS, + .init = km_init_kmeans, +}}; + +#define NUM_TYPES (sizeof(TYPES) / sizeof(TYPES[0])) + +km_init_type_t +km_init_get_type( + const char * const s +) { + // find init method + for (size_t i = 0; i < NUM_TYPES; i++) { + if (!strcmp(s, TYPES[i].name)) { + // return type + return TYPES[i].type; + } + } + + // return failure + return KM_INIT_TYPE_LAST; +} + +static km_init_fn_t +get_init_func( + const km_init_type_t type +) { + // find init method + for (size_t i = 0; i < NUM_TYPES; i++) { + if (type == TYPES[i].type) { + // return init func + return TYPES[i].init; + } + } + + // return failure + return NULL; +} + +bool +km_init( + km_set_t * const cs, + const km_init_type_t init_type, + const size_t num_clusters, + const km_set_t * const set, + km_rand_t * const rs +) { + km_init_fn_t init_fn = get_init_func(init_type); + if (!init_fn) { + die("unknown cluster init method"); + } + + // call init, return result + return init_fn(cs, num_clusters, set, rs); +} diff --git a/src/km-load.c b/src/km-load.c new file mode 100644 index 0000000..3615aa9 --- /dev/null +++ b/src/km-load.c @@ -0,0 +1,112 @@ +#include // bool +#include // fscanf() +#include // strerror() +#include // errno +#include "util.h" +#include "km.h" + +#define FAIL(...) do { \ + if (cbs && cbs->on_error) { \ + char buf[1024]; \ + snprintf(buf, sizeof(buf), __VA_ARGS__); \ + cbs->on_error(buf, cb_data); \ + } \ + return false; \ +} while (0) + +_Bool +km_load( + FILE * const fh, + const km_load_cbs_t * const cbs, + void * const cb_data +) { + // read shape + km_shape_t shape = { 0, 0 }; + if (fscanf(fh, "%zu %zu", &(shape.num_floats), &(shape.num_ints)) != 2) { + FAIL("shape fscanf() failed: %s", strerror(errno)); + } + + if (cbs && cbs->on_shape) { + // emit shape + if (!cbs->on_shape(&shape, cb_data)) { + // return failure + return false; + } + } + + // alloc floats buffer + float *floats = NULL; + if (shape.num_floats > 0) { + floats = malloc(sizeof(float) * shape.num_floats); + if (!floats) { + FAIL("floats malloc() failed: %s", strerror(errno)); + } + } + + // alloc ints buffer + int *ints = NULL; + if (shape.num_ints > 0) { + ints = malloc(sizeof(int) * shape.num_ints); + if (!ints) { + FAIL("ints malloc() failed: %s", strerror(errno)); + } + } + + for (size_t row = 0; !feof(fh); row++) { + for (size_t i = 0; i < shape.num_floats; i++) { + if (fscanf(fh, " %f ", floats + i) != 1) { + FAIL("[%zu, %zu] float fscanf() failed: %s", row, i, strerror(errno)); + } + } + + // read ints + for (size_t i = 0; i < shape.num_ints; i++) { + if (fscanf(fh, " %d ", ints + i) != 1) { + FAIL("[%zu, %zu] int fscanf() failed: %s", row, i, strerror(errno)); + } + } + + if (cbs && cbs->on_row) { + // emit row + if (!cbs->on_row(floats, ints, cb_data)) { + // return failure + return false; + } + } + } + + if (shape.num_floats > 0) { + // free float buffer + free(floats); + } + + if (shape.num_ints > 0) { + // free int buffer + free(ints); + } + + // return success + return true; +} + +_Bool +km_load_path( + const char * const path, + const km_load_cbs_t * const cbs, + void * const cb_data +) { + // open file + FILE *fh = fopen(path, "rb"); + if (!fh) { + FAIL("fopen(\"%s\") failed: %s", path, strerror(errno)); + } + + // load file, get result + const bool r = km_load(fh, cbs, cb_data); + + // close file + fclose(fh); + + // return result + return r; +} diff --git a/src/km-rand-erand48.c b/src/km-rand-erand48.c new file mode 100644 index 0000000..de84b93 --- /dev/null +++ b/src/km-rand-erand48.c @@ -0,0 +1,89 @@ +#include // bool +#define _DEFAULT_SOURCE +#include // drand48_r() +#include // fopen() +#include "util.h" +#include "km.h" + +// get N get_floats +static bool +on_get_floats( + km_rand_t * const rs, + const size_t num_vals, + float * const vals +) { + unsigned short * const data = rs->data; + + // generate results + for (size_t i = 0; i < num_vals; i++) { + vals[i] = erand48(data); + } + + // return success + return true; +} + +// fill sizes callback for system random source +static bool +on_get_sizes( + km_rand_t * const rs, + const size_t num_vals, + size_t * const vals +) { + unsigned short * const data = rs->data; + + // generate results + for (size_t i = 0; i < num_vals; i++) { + vals[i] = nrand48(data); + } + + // return success + return true; +} + +static void +on_fini( + km_rand_t * const rs +) { + unsigned short * data = rs->data; + + // free rng data + free(data); + + // clear pointer + rs->data = NULL; +} + +// path random source callbacks +static const km_rand_cbs_t +ERAND48_RAND_CBS = { + .get_floats = on_get_floats, + .get_sizes = on_get_sizes, + .fini = on_fini, +}; + +// init system random source (uses system rand()) +bool +km_rand_init_erand48( + km_rand_t * const rs, + const uint64_t seed +) { + // alloc rng data + unsigned short * const data = malloc(3 * sizeof(unsigned short)); + if (!data) { + // return failure + return false; + } + + // populate seed + data[0] = seed & 0xffff; + data[1] = (seed >> 16) & 0xffff; + data[2] = (seed >> 32) & 0xffff; + + // populate context + rs->cbs = &ERAND48_RAND_CBS; + rs->data = data; + + // return success + return true; +} diff --git a/src/km-rand-libc.c b/src/km-rand-libc.c new file mode 100644 index 0000000..da1e4b0 --- /dev/null +++ b/src/km-rand-libc.c @@ -0,0 +1,56 @@ +#include // bool +#include "util.h" +#include "km.h" + +// libc random source get_floats +static bool +on_get_floats( + km_rand_t * const rs, + const size_t num_vals, + float * const vals +) { + UNUSED(rs); + + // generate results + for (size_t i = 0; i < num_vals; i++) { + vals[i] = 1.0 * rand() / RAND_MAX; + } + + // return success + return true; +} + +// fill sizes callback for system random source +static bool +on_get_sizes( + km_rand_t * const rs, + const size_t num_vals, + size_t * const vals +) { + UNUSED(rs); + + // generate results + for (size_t i = 0; i < num_vals; i++) { + vals[i] = rand(); + } + + // return success + return true; +} + +// system random source callbacks +static const km_rand_cbs_t +LIBC_RAND_CBS = { + .get_floats = on_get_floats, + .get_sizes = on_get_sizes, + .fini = NULL, +}; + +// init system random source (uses system rand()) +void +km_rand_init_libc( + km_rand_t * const rs +) { + rs->cbs = &LIBC_RAND_CBS; + rs->data = NULL; +} diff --git a/src/km-rand-path.c b/src/km-rand-path.c new file mode 100644 index 0000000..6597038 --- /dev/null +++ b/src/km-rand-path.c @@ -0,0 +1,152 @@ +#include // bool +#include // memset() +#include // fopen() +#include "util.h" +#include "km.h" + +#define MAX_BUF_ITEMS (1 << 20) + +typedef struct { + FILE *fh; + size_t ofs; + uint64_t buf[MAX_BUF_ITEMS]; +} ctx_t; + +static inline bool +refill( + ctx_t * const ctx +) { + ctx->ofs = 0; + // D("reading %zu bytes", sizeof(ctx->buf)); + return !!fread(ctx->buf, sizeof(ctx->buf), 1, ctx->fh); +} + +static bool +next( + ctx_t * const ctx, + uint64_t * const r +) { + if ((ctx->ofs >= MAX_BUF_ITEMS - 1) && !refill(ctx)) { + // return failure + return false; + } + + // get next value, increment offset + *r = ctx->buf[ctx->ofs]; + ctx->ofs++; + + // return success + return true; +} + +// get N get_floats +static bool +on_get_floats( + km_rand_t * const rs, + const size_t num_vals, + float * const vals +) { + ctx_t * const ctx = rs->data; + + // generate results + for (size_t i = 0; i < num_vals; i++) { + uint64_t v; + + if (!next(ctx, &v)) { + // return failure + return false; + } + + // add value + vals[i] = 1.0 * v / RAND_MAX; + } + + // return success + return true; +} + +// fill sizes callback for system random source +static bool +on_get_sizes( + km_rand_t * const rs, + const size_t num_vals, + size_t * const vals +) { + ctx_t * const ctx = rs->data; + + // generate results + for (size_t i = 0; i < num_vals; i++) { + uint64_t v; + + if (!next(ctx, &v)) { + // return failure + return false; + } + + // add value + vals[i] = v; + } + + // return success + return true; +} + +static void +on_fini( + km_rand_t * const rs +) { + ctx_t * const ctx = rs->data; + + // close file + fclose(ctx->fh); + + // free context + free(ctx); + + rs->data = NULL; +} + +// path random source callbacks +static const km_rand_cbs_t +PATH_RAND_CBS = { + .get_floats = on_get_floats, + .get_sizes = on_get_sizes, + .fini = on_fini, +}; + +// init system random source (uses system rand()) +bool +km_rand_init_path( + km_rand_t * const rs, + const char * const path +) { + // open file + FILE *fh = fopen(path, "rb"); + if (!fh) { + // return failure + return false; + } + + // alloc context + ctx_t * const ctx = malloc(sizeof(ctx_t)); + if (!ctx) { + // return failure + return false; + } + + // populate context + ctx->fh = fh; + + // fill buffer + if (!refill(ctx)) { + // return failure + return false; + } + + // populate random source + rs->cbs = &PATH_RAND_CBS; + rs->data = ctx; + + // return success + return true; +} diff --git a/src/km-rand.c b/src/km-rand.c new file mode 100644 index 0000000..e777069 --- /dev/null +++ b/src/km-rand.c @@ -0,0 +1,33 @@ +#include +#include "util.h" +#include "km.h" + +// fill buffer with N random floats +bool +km_rand_get_floats( + km_rand_t * const rs, + const size_t num_floats, + float * const floats +) { + return rs->cbs->get_floats(rs, num_floats, floats); +} + +// fill buffer with N random size_ts +bool +km_rand_get_sizes( + km_rand_t * const rs, + const size_t num_sizes, + size_t * const sizes +) { + return rs->cbs->get_sizes(rs, num_sizes, sizes); +} + +// finalize random source +void +km_rand_fini( + km_rand_t * const rs +) { + if (rs->cbs->fini) { + rs->cbs->fini(rs); + } +} diff --git a/src/km-set.c b/src/km-set.c new file mode 100644 index 0000000..69ae6c0 --- /dev/null +++ b/src/km-set.c @@ -0,0 +1,267 @@ +#include // bool +#include // size_t +#include // memcpy() +#include // rand() +#include // errno +#include "util.h" +#include "km.h" + +#define MIN_ROWS (4096 / sizeof(float)) + +// grow data set +static bool +km_set_grow( + km_set_t * const set, + const size_t capacity +) { + float *floats = NULL; + const size_t floats_size = sizeof(float) * set->shape.num_floats * capacity; + // fprintf(stderr, "floats_size = %zu\n", floats_size); + if (floats_size > 0) { + // alloc floats + floats = realloc(set->floats, floats_size); + if (!floats) { + // return failure + return false; + } + } + + int *ints = NULL; + const size_t ints_size = sizeof(int) * set->shape.num_ints * capacity; + // fprintf(stderr, "ints_size = %zu\n", ints_size); + if (ints_size > 0) { + // alloc ints + ints = realloc(set->ints, ints_size); + if (!ints) { + // return failure + return false; + } + } + + // update set + set->floats = floats; + set->ints = ints; + set->capacity = capacity; + + // return success + return true; +} + +// init data set with shape and initial size +bool +km_set_init( + km_set_t * const set, + const km_shape_t * const shape, + const size_t row_capacity +) { + // alloc bounds + float * const bounds = malloc(2 * sizeof(float) * shape->num_floats); + if (!bounds) { + // return error + return false; + } + + set->state = KM_SET_STATE_INIT; + set->floats = NULL; + set->ints = NULL; + set->shape = *shape; + set->num_rows = 0; + set->capacity = 0; + set->bounds = bounds; + + return km_set_grow(set, MAX(MIN_ROWS, row_capacity + 1)); +} + +// finalize data set +void +km_set_fini(km_set_t * const set) { + if (set->state == KM_SET_STATE_FINI) { + return; + } + + if (set->bounds) { + // free bounds + free(set->bounds); + set->bounds = NULL; + } + + if (set->floats) { + // free floats + free(set->floats); + set->floats = NULL; + } + + if (set->ints) { + // free ints + free(set->ints); + set->ints = NULL; + } + + // shrink capacity + set->capacity = 0; + + // set state + set->state = KM_SET_STATE_FINI; +} + +// append rows to data set, growing set if necessary +bool +km_set_push( + km_set_t * const set, + const size_t num_rows, + const float * const floats, + const int * const ints +) { + // check state + if (set->state != KM_SET_STATE_INIT) { + // return failure + return false; + } + + const size_t capacity_needed = set->num_rows + num_rows; + // FIXME: potential overflow here + if (capacity_needed >= set->capacity) { + // crappy growth algorithm + const size_t new_capacity = 2 * capacity_needed + 1; + + // resize storage + if (!km_set_grow(set, MAX(MIN_ROWS, new_capacity))) { + return false; + } + } + + // copy floats + const size_t num_floats = set->shape.num_floats; + if (num_floats > 0) { + float * const dst = set->floats + num_floats * set->num_rows; + const size_t stride = sizeof(float) * num_floats; + + // copy floats + memcpy(dst, floats, stride * num_rows); + + if (!set->num_rows) { + // there were no rows, so populate bounds with first row + memcpy(set->bounds, floats, stride); + memcpy(set->bounds + num_floats, floats, stride); + } + + for (size_t i = 0; i < num_rows; i++) { + for (size_t j = 0; j < num_floats; j++) { + const float val = floats[i * num_floats + j]; + + if (val < set->bounds[j]) { + // update min bound + set->bounds[j] = val; + } + + if (val > set->bounds[num_floats + j]) { + // update max bound + set->bounds[num_floats + j] = val; + } + } + } + } + + // copy ints + const size_t num_ints = set->shape.num_ints; + if (num_ints > 0) { + int * const dst = set->ints + num_ints * set->num_rows; + const size_t stride = sizeof(int) * num_ints; + + // copy ints + memcpy(dst, ints, stride * num_rows); + } + + // increment row count + set->num_rows += num_rows; + + // return success + return true; +} + +bool +km_set_copy( + km_set_t * const dst, + const km_set_t * const src +) { + if (src->state != KM_SET_STATE_INIT && src->state != KM_SET_STATE_NORMALIZED) { + // return failure + D("invalid state"); + return false; + } + + // init dst set + if (!km_set_init(dst, &(src->shape), src->num_rows)) { + // return failure + D("km_set_init()"); + return false; + } + + // copy floats + if (src->shape.num_floats > 0) { + const size_t stride = sizeof(float) * src->shape.num_floats; + + // copy floats + memcpy(dst->floats, src->floats, stride * src->num_rows); + + // copy bounds + memcpy(dst->bounds, src->bounds, 2 * stride); + } + + // copy ints + const size_t num_ints = src->shape.num_ints; + if (num_ints > 0) { + const size_t stride = sizeof(int) * num_ints; + + // copy ints + memcpy(dst->ints, src->ints, stride * src->num_rows); + } + + // increment row count + dst->num_rows = src->num_rows; + + // return success + return true; +} + +bool +km_set_normalize( + km_set_t * const set +) { + const size_t num_floats = set->shape.num_floats; + + // check set state + if (set->state != KM_SET_STATE_INIT) { + // return failure + return false; + } + + // normalize values + for (size_t i = 0; i < set->num_rows; i++) { + for (size_t j = 0; j < num_floats; j++) { + const size_t ofs = i * num_floats + j; + const float val = set->floats[ofs], + min = set->bounds[j], + max = set->bounds[num_floats + j]; + + // normalize and write value + set->floats[ofs] = (val - min) / (max - min); + } + } + + // set state + set->state = KM_SET_STATE_NORMALIZED; + + // return success + return true; +} + +// get row from data set +float * +km_set_get_row( + const km_set_t * const set, + const size_t i +) { + const size_t num_floats = set->shape.num_floats; + return set->floats + i * num_floats; +} diff --git a/src/km-solve.c b/src/km-solve.c new file mode 100644 index 0000000..f579a41 --- /dev/null +++ b/src/km-solve.c @@ -0,0 +1,200 @@ +#include // bool +#include // memset() +#include // FLT_MAX +#include // sqrt() +#include "util.h" +#include "km.h" + +#define MIN_CLUSTER_DISTANCE 0.00001 + +// alloc and initialize row map +static km_row_map_t * +km_row_map_init( + const size_t num_rows +) { + // alloc row map + km_row_map_t * const row_map = malloc(sizeof(km_row_map_t) * num_rows); + + // check for error + if (!row_map) { + // return failure + return false; + } + + // init row map + for (size_t i = 0; i < num_rows; i++) { + // setting distances to maximum + row_map[i].d2 = FLT_MAX; + row_map[i].d2_near = FLT_MAX; + } + + // return row map + return row_map; +} + +static void +km_row_map_fini( + km_row_map_t * const row_map +) { + free(row_map); +} + +// use k-means to iteratively update the cluster centroids until there +// are no more changes to the centroids +bool +km_solve( + km_set_t * const cs, + const km_set_t * const set, + const km_solve_cbs_t * const cbs, + void * const cb_data +) { + const size_t num_clusters = cs->num_rows, + num_floats = set->shape.num_floats; + + // row map: row => distance, cluster ID + km_row_map_t *row_map = km_row_map_init(set->num_rows); + if (!row_map) { + // return failure + return false; + } + + // calculate clusters by doing the following: + // * walk all clusters and all rows + // * if we find a closer cluster, move row to cluster + // * if there were changes to any cluster, then calculate a new + // centroid for each cluster by averaging the cluster rows + // * repeat until there are no more changes + bool changed = false; + do { + // no changes yet + changed = false; + + for (size_t i = 0; i < num_clusters; i++) { + // get the floats for this cluster + const float * const floats = km_set_get_row(cs, i); + + for (size_t j = 0; j < set->num_rows; j++) { + // get row values + const float * const vals = km_set_get_row(set, j); + + // calculate the distance squared between row and cluster + const float d2 = distance_squared(num_floats, floats, vals); + + if (d2 < row_map[j].d2) { + // row is closer to this cluster, update row map + row_map[j].d2 = d2; + row_map[j].cluster = i; + + // flag change + changed = true; + } + + if ((row_map[j].cluster != i) && (d2 < row_map[j].d2_near)) { + row_map[j].d2_near = d2; + + // flag change + changed = true; + } + } + } + + if (changed) { + // if there were changes, then we need to calculate the new + // cluster centers + + // calculate new center + for (size_t i = 0; i < num_clusters; i++) { + size_t num_rows = 0; + float * const floats = km_set_get_row(cs, i); + memset(floats, 0, sizeof(float) * num_floats); + + for (size_t j = 0; j < set->num_rows; j++) { + const float * const row_floats = km_set_get_row(set, j); + + if (row_map[j].cluster == i) { + // calculate numerator for average + for (size_t k = 0; k < num_floats; k++) { + floats[k] += row_floats[k]; + } + + // increment denominator for average + num_rows++; + } + } + + // save number of rows in this cluster + cs->ints[i] = num_rows; + + if (num_rows > 0) { + for (size_t k = 0; k < num_floats; k++) { + // divide by denominator to get average + floats[k] /= num_rows; + } + } + } + } + + if (cbs && cbs->on_step) { + // pass clusters and row map to step callback + cbs->on_step(cs, row_map, cb_data); + } + } while (changed); + + if (cbs && cbs->on_stats) { + float sum = 0, + silouette = 0, + mean_dists[num_clusters], + mean_nears[num_clusters]; + + memset(mean_dists, 0, sizeof(mean_dists)); + memset(mean_nears, 0, sizeof(mean_nears)); + + // calculate sum of distances across all clusters + for (size_t i = 0; i < set->num_rows; i++) { + sum += row_map[i].d2; + } + + // calculate mean numerators and silouette + for (size_t i = 0; i < set->num_rows; i++) { + // distance squared (d2) to center of this cluster + mean_dists[row_map[i].cluster] += row_map[i].d2; + + // distance squared (d2) to center of nearest cluster + mean_nears[row_map[i].cluster] += row_map[i].d2_near; + + // calculate silouette denominator + // (https://en.wikipedia.org/wiki/Silhouette_%28clustering%29) + const float delta = row_map[i].d2_near - row_map[i].d2; + if (fabsf(delta) > MIN_CLUSTER_DISTANCE) { + silouette += delta / MAX(row_map[i].d2, row_map[i].d2_near); + } + } + + // finalize means (divide by row count) + for (size_t i = 0; i < num_clusters; i++) { + mean_dists[i] = (cs->ints[i]) ? (sqrt(mean_dists[i]) / cs->ints[i]) : 0; + mean_nears[i] = (cs->ints[i]) ? (sqrt(mean_nears[i]) / cs->ints[i]) : 0; + } + + // finalize silouette + silouette /= set->num_rows; + + // build stats + const km_solve_stats_t stats = { + .sum = sum, + .silouette = silouette, + .mean_dists = mean_dists, + .mean_nears = mean_nears, + .num_clusters = num_clusters, + }; + + // emit means + cbs->on_stats(cs, &stats, cb_data); + } + + // free row_map + km_row_map_fini(row_map); + + // return success + return true; +} diff --git a/src/km.h b/src/km.h new file mode 100644 index 0000000..5a34261 --- /dev/null +++ b/src/km.h @@ -0,0 +1,230 @@ +#ifndef KM_H +#define KM_H + +#include // size_t +#include // uint8_t +#include // FILE + +// forward typedef for callbacks +typedef struct km_rand_t_ km_rand_t; + +// random number source callbacks +typedef struct { + _Bool (*get_floats)(km_rand_t * const, const size_t, float * const); + _Bool (*get_sizes)(km_rand_t * const, const size_t, size_t * const); + void (*fini)(km_rand_t * const); +} km_rand_cbs_t; + +// random number source +struct km_rand_t_ { + const km_rand_cbs_t *cbs; + void *data; +}; + +// fill buffer with N random floats +_Bool km_rand_get_floats(km_rand_t * const, const size_t, float * const); + +// fill buffer with N random size_ts +_Bool km_rand_get_sizes(km_rand_t * const, const size_t, size_t * const); + +// finalize random source +void km_rand_fini(km_rand_t * const); + +// init libc rng (use libc rand()) +void km_rand_init_libc(km_rand_t *); + +// init file source rng (e.g. "/dev/urandom") +_Bool km_rand_init_path(km_rand_t * const, const char *); + +// init erand48 rng (use erand48()) +_Bool km_rand_init_erand48(km_rand_t * const, const uint64_t); + +// shape of data set +typedef struct { + size_t num_floats, + num_ints; +} km_shape_t; + +typedef enum { + KM_SET_STATE_INIT, + KM_SET_STATE_NORMALIZED, + KM_SET_STATE_FINI, + KM_SET_STATE_LAST, +} km_set_state_t; + +// data set +typedef struct { + km_set_state_t state; + + float *floats, + *bounds; + int *ints; + + km_shape_t shape; + + size_t num_rows, + capacity; +} km_set_t; + +// init data set with shape and initial size +_Bool km_set_init(km_set_t * const, const km_shape_t *, const size_t); +// finalize data set +void km_set_fini(km_set_t * const); + +// append rows to data set, growing set if necessary +_Bool km_set_push(km_set_t *, const size_t, const float *, const int *); + +// deep copy data set +_Bool km_set_copy(km_set_t * const, const km_set_t * const); + +// normalize data set +_Bool km_set_normalize(km_set_t * const); + +// get pointer to data set row +float *km_set_get_row(const km_set_t *, const size_t); + +typedef struct { + float d2, + d2_near; + size_t cluster; +} km_row_map_t; + +typedef struct { + const float sum; + + const float silouette; + + const float *mean_dists; + const float *mean_nears; + + const size_t num_clusters; +} km_solve_stats_t; + +typedef struct { + void (*on_step)(const km_set_t *, const km_row_map_t *, void *); + void (*on_stats)(const km_set_t *, const km_solve_stats_t *, void *); +} km_solve_cbs_t; + +// use k-means to iteratively update the cluster centroids until there +// are no more changes to the centroids +_Bool +km_solve( + km_set_t * const, + const km_set_t * const, + const km_solve_cbs_t * const, + void * const +); + +typedef struct { + const km_set_t * const cluster_set; + + float distance_sum, + silouette, + mean_cluster_size; + + const size_t num_clusters, + num_empty_clusters; +} km_find_data_t; + +typedef _Bool (*km_find_init_cb_t)( + km_set_t * const cluster_set, + const size_t num_clusters, + const km_set_t * const data_set, + void *cb_data +); + +typedef _Bool (*km_find_fini_cb_t)( + km_set_t * const cluster_set, + void *cb_data +); + +typedef _Bool (*km_find_best_cb_t)( + const float score, + const km_set_t * const cluster_set, + void *cb_data +); + +typedef void (*km_find_data_cb_t)( + const km_find_data_t * const, + void * +); + +typedef struct { + // maximum number of clusters + size_t max_clusters; + + // maximum number of tests per cluster + size_t num_tests; + + km_find_init_cb_t on_init; + km_find_fini_cb_t on_fini; + km_find_data_cb_t on_data; + km_find_best_cb_t on_best; +} km_find_cbs_t; + +// repeatedly test different cluster sizes and report results +_Bool +km_find( + const km_set_t * const set, + const km_find_cbs_t * const cbs, + void * const cb_data +); + +// render set to rgb buffer +void +km_set_draw( + const km_set_t * const set, + uint8_t * const rgb, + const size_t width, + const size_t height, + const int dot_size, + const uint32_t color +); + +// load callbacks +typedef struct { + _Bool (*on_shape)(const km_shape_t * const, void *); + _Bool (*on_row)(const float *, const int *, void *); + void (*on_error)(const char * const, void *); +} km_load_cbs_t; + +// load set from file handle +_Bool +km_load( + FILE * const fh, + const km_load_cbs_t * const cbs, + void * const cb_data +); + +// load set from path +_Bool km_load_path( + const char * const path, + const km_load_cbs_t * const cbs, + void * const cb_data +); + +typedef enum { + // rand: init cluster set by picking random centroids + KM_INIT_TYPE_RAND, + + // forgy: init cluster set by picking random rows from set + KM_INIT_TYPE_FORGY, + + // kmeans++: ... + KM_INIT_TYPE_KMEANS, + + KM_INIT_TYPE_LAST, +} km_init_type_t; + +km_init_type_t km_init_get_type(const char * const); + +// init cluster set using given method +_Bool km_init( + km_set_t * const, + const km_init_type_t, + const size_t, + const km_set_t * const, + km_rand_t * const +); + +#endif /* KM_H */ diff --git a/src/main.c b/src/main.c new file mode 100644 index 0000000..55056b7 --- /dev/null +++ b/src/main.c @@ -0,0 +1,398 @@ +#include // bool +#include // fprintf() +#include // EXIT_{FAILURE,SUCCESS} +#include // exit() +#include // memset() +#include // fabsf() + +#define STB_IMAGE_WRITE_IMPLEMENTATION +#define STB_ONLY_PNG +#include "stb_image_write.h" + +#include "util.h" +#include "km.h" + +#define MAX_CLUSTERS 10 +#define NUM_TESTS 300 +#define MAX_BEST 10 + +#define IM_WIDTH 128 +#define IM_HEIGHT 128 +#define IM_STRIDE (3 * IM_WIDTH) + +typedef struct { + float score; + km_set_t set; +} best_item_t; + +typedef struct { + // cluster initialization method + km_init_type_t init_type; + + // random number source + km_rand_t rs; + + struct { + float distance, + silouette, + cluster_size; + size_t num_empty; + } rows[MAX_CLUSTERS - 2]; + + // best clusters + best_item_t best[MAX_BEST]; + size_t num_best; +} ctx_t; + +static int +best_score_cmp( + const void * const ap, + const void * const bp +) { + const best_item_t * const a = ap; + const best_item_t * const b = bp; + + return (a->score > b->score) ? 1 : -1; +} + +static void +ctx_best_sort( + ctx_t * const ctx +) { + // sort best sets by ascending score (worst to best) + qsort( + ctx->best, + MIN(ctx->num_best, MAX_BEST), + sizeof(best_item_t), + best_score_cmp + ); +} + +static void +ctx_best_walk( + const ctx_t * const ctx, + void (*on_best)(const km_set_t * const, const size_t, const float, void *), + void * const cb_data +) { + if (on_best) { + // walk best sets and emit each one + for (size_t i = 0; i < MIN(ctx->num_best, MAX_BEST); i++) { + on_best(&(ctx->best[i].set), i, ctx->best[i].score, cb_data); + } + } +} + +static bool +load_on_shape( + const km_shape_t * const shape, + void * const cb_data +) { + km_set_t * const set = cb_data; + + // D("shape: %zu floats, %zu ints", shape->num_floats, shape->num_ints); + + // init set + if (!km_set_init(set, shape, 100)) { + die("km_set_init() failed"); + } + + // return success + return true; +} + +static bool +load_on_row( + const float * const floats, + const int * const ints, + void * const cb_data +) { + km_set_t * const set = cb_data; + + // push row + if (!km_set_push(set, 1, floats, ints)) { + die("km_set_push_rows() failed"); + } + + // return success + return true; +} + +static void +load_on_error( + const char * const err, + void * const cb_data +) { + UNUSED(cb_data); + die("load failed: %s", err); +} + +static const km_load_cbs_t +LOAD_CBS = { + .on_shape = load_on_shape, + .on_row = load_on_row, + .on_error = load_on_error, +}; + +static bool +find_on_init( + km_set_t * const cs, + const size_t num_clusters, + const km_set_t * const set, + void *cb_data +) { + ctx_t * const ctx = cb_data; + return km_init(cs, ctx->init_type, num_clusters, set, &(ctx->rs)); +} + +static bool +find_on_fini( + km_set_t * const cs, + void *cb_data +) { + UNUSED(cb_data); + km_set_fini(cs); + return true; +} + +static void +find_on_data( + const km_find_data_t * const data, + void *cb_data +) { + ctx_t * const ctx = cb_data; + const size_t ofs = data->num_clusters - 2; + + ctx->rows[ofs].distance += data->distance_sum; + ctx->rows[ofs].silouette += data->silouette; + ctx->rows[ofs].cluster_size += data->mean_cluster_size; + ctx->rows[ofs].num_empty += data->num_empty_clusters; +} + +static bool +find_on_best( + const float score, + const km_set_t * const cs, + void *cb_data +) { + ctx_t * const ctx = cb_data; + + D("new best: score = %0.3f, num_clusters = %zu", score, cs->num_rows); + + // get pointer to destination set + // (note: data->best is a ring buffer) + const size_t ofs = ctx->num_best % MAX_BEST; + ctx->best[ofs].score = score; + km_set_t * const dst = &(ctx->best[ofs].set); + + if (ctx->num_best >= MAX_BEST) { + // finalize old best data set + km_set_fini(dst); + } + + // copy data set to best ring buffer + if (!km_set_copy(dst, cs)) { + die("km_set_copy()"); + } + + // increment best count + ctx->num_best++; + + // return success + return true; +} + +// init find config +static const km_find_cbs_t +FIND_CBS = { + .max_clusters = MAX_CLUSTERS, + .num_tests = NUM_TESTS, + + .on_init = find_on_init, + .on_fini = find_on_fini, + .on_data = find_on_data, + .on_best = find_on_best, +}; + +static void +ctx_csv_print_row( + const ctx_t * const ctx, + FILE * const fh, + const size_t i +) { + const size_t num_clusters = i + 2; + const float mean_distance = ctx->rows[i].distance / NUM_TESTS, + mean_cluster_size = ctx->rows[i].cluster_size / NUM_TESTS, + mean_empty = 1.0 * ctx->rows[i].num_empty / NUM_TESTS, + score = ctx->rows[i].silouette / NUM_TESTS; + + // print result + fprintf(fh, "%zu,%0.3f,%0.3f,%0.3f,%0.3f\n", + num_clusters, + score, + mean_distance, + mean_cluster_size, + mean_empty + ); +} + +static void +ctx_csv_print( + const ctx_t * const ctx, + FILE * const fh +) { + // print headers + fprintf(fh, "#,score,distance,size,empty\n"); + + // print rows + for (size_t i = 0; i < MAX_CLUSTERS - 2; i++) { + ctx_csv_print_row(ctx, fh, i); + } +} + +// static image data buffer +static uint8_t im_data[3 * IM_WIDTH * IM_HEIGHT]; + +static void +save_on_best( + const km_set_t * const set, + const size_t rank, + const float score, + void * const cb_data +) { + const ctx_t * const ctx = cb_data; + UNUSED(score); + + // convert rank to channel brightness + const uint8_t ch = 0x33 + (0xff - 0x33) * (1.0 * rank + 1) / (MAX_BEST); + const uint8_t shift = (rank == MIN(ctx->num_best, MAX_BEST) - 1) ? 8 : 16; + const uint32_t color = (ch & 0xff) << shift; + // const uint32_t color = 0xff0000; + // D("rank = %zu, score = %0.3f, size = %zu, color = %06x", rank, score, set->num_rows, color); + + // draw clusters + km_set_draw(set, im_data, IM_WIDTH, IM_HEIGHT, 3, color); +} + +static void +ctx_save_png( + const ctx_t * const ctx, + const char * const png_path, + const km_set_t * const set +) { + // clear image data to white + memset(im_data, 0xff, sizeof(im_data)); + + // draw data points + km_set_draw(set, im_data, IM_WIDTH, IM_HEIGHT, 1, 0x000000); + if (!stbi_write_png(png_path, IM_WIDTH, IM_HEIGHT, 3, im_data, IM_STRIDE)) { + die("stbi_write_png(\"%s\")", png_path); + } + + // draw best cluster points + ctx_best_walk(ctx, save_on_best, (void*) ctx); + + // save png + if (!stbi_write_png(png_path, IM_WIDTH, IM_HEIGHT, 3, im_data, IM_STRIDE)) { + die("stbi_write_png(\"%s\")", png_path); + } +} + +static void +ctx_best_print_on_best( + const km_set_t * const set, + const size_t rank, + const float score, + void * const cb_data +) { + FILE * const fh = cb_data; + + fprintf(fh, "rank = %zu, score = %0.3f, num_clusters = %zu: [\n", rank, score, set->num_rows); + for (size_t i = 0; i < set->num_rows; i++) { + const float * const vals = km_set_get_row(set, i); + fprintf(fh, " ["); + + for (size_t j = 0; j < set->shape.num_floats; j++) { + fprintf(fh, "%s%0.3f", (j > 0) ? ", " : "", vals[j]); + } + + fprintf(fh, "], (%d rows)\n", set->ints[i]); + } + fprintf(fh, "]\n"); +} + +static void +ctx_best_print( + const ctx_t * const ctx, + FILE * const fh +) { + ctx_best_walk(ctx, ctx_best_print_on_best, fh); +} + +static const char USAGE_FORMAT[] = + "Usage: %s [init] [data_path] \n" + "\n" + "Arguments:\n" + "* init: Cluster init method (one of \"rand\" or \"set\").\n" + "* data_path: Path to input data file.\n" + "* png_path: Path to output file (optional).\n" + ""; + +int main(int argc, char *argv[]) { + // check command-line arguments + if (argc < 3) { + fprintf(stderr, USAGE_FORMAT, argv[0]); + return EXIT_FAILURE; + } + + // get command-line arguments + const char * const init_type_name = argv[1]; + const char * const data_path = argv[2]; + const char * const png_path = (argc > 3) ? argv[3] : NULL; + + // init random seed + srand(getpid()); + + // init context + ctx_t ctx; + memset(&ctx, 0, sizeof(ctx_t)); + ctx.init_type = km_init_get_type(init_type_name); + + // init ctx rng + km_rand_init_erand48(&(ctx.rs), rand()); + + // init data set + km_set_t set; + if (!km_load_path(data_path, &LOAD_CBS, &set)) { + die("km_load_path(\"%s\") failed", data_path); + } + + // normalize data set + if (!km_set_normalize(&set)) { + die("km_set_normalize() failed"); + } + + // find best solutions + if (!km_find(&set, &FIND_CBS, &ctx)) { + die("km_find()"); + } + + // print csv + ctx_csv_print(&ctx, stdout); + + // sort best results from lowest to highest + ctx_best_sort(&ctx); + + // print best + ctx_best_print(&ctx, stdout); + + if (png_path) { + // save png of normalized data set and best clusters + ctx_save_png(&ctx, png_path, &set); + } + + // finalize data set + km_set_fini(&set); + + // return success + return 0; +} diff --git a/src/stb_image_write.h b/src/stb_image_write.h new file mode 100644 index 0000000..c05e958 --- /dev/null +++ b/src/stb_image_write.h @@ -0,0 +1,1568 @@ +/* stb_image_write - v1.09 - public domain - http://nothings.org/stb/stb_image_write.h + writes out PNG/BMP/TGA/JPEG/HDR images to C stdio - Sean Barrett 2010-2015 + no warranty implied; use at your own risk + + Before #including, + + #define STB_IMAGE_WRITE_IMPLEMENTATION + + in the file that you want to have the implementation. + + Will probably not work correctly with strict-aliasing optimizations. + + If using a modern Microsoft Compiler, non-safe versions of CRT calls may cause + compilation warnings or even errors. To avoid this, also before #including, + + #define STBI_MSC_SECURE_CRT + +ABOUT: + + This header file is a library for writing images to C stdio. It could be + adapted to write to memory or a general streaming interface; let me know. + + The PNG output is not optimal; it is 20-50% larger than the file + written by a decent optimizing implementation; though providing a custom + zlib compress function (see STBIW_ZLIB_COMPRESS) can mitigate that. + This library is designed for source code compactness and simplicity, + not optimal image file size or run-time performance. + +BUILDING: + + You can #define STBIW_ASSERT(x) before the #include to avoid using assert.h. + You can #define STBIW_MALLOC(), STBIW_REALLOC(), and STBIW_FREE() to replace + malloc,realloc,free. + You can #define STBIW_MEMMOVE() to replace memmove() + You can #define STBIW_ZLIB_COMPRESS to use a custom zlib-style compress function + for PNG compression (instead of the builtin one), it must have the following signature: + unsigned char * my_compress(unsigned char *data, int data_len, int *out_len, int quality); + The returned data will be freed with STBIW_FREE() (free() by default), + so it must be heap allocated with STBIW_MALLOC() (malloc() by default), + +USAGE: + + There are five functions, one for each image file format: + + int stbi_write_png(char const *filename, int w, int h, int comp, const void *data, int stride_in_bytes); + int stbi_write_bmp(char const *filename, int w, int h, int comp, const void *data); + int stbi_write_tga(char const *filename, int w, int h, int comp, const void *data); + int stbi_write_jpg(char const *filename, int w, int h, int comp, const void *data, int quality); + int stbi_write_hdr(char const *filename, int w, int h, int comp, const float *data); + + void stbi_flip_vertically_on_write(int flag); // flag is non-zero to flip data vertically + + There are also five equivalent functions that use an arbitrary write function. You are + expected to open/close your file-equivalent before and after calling these: + + int stbi_write_png_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data, int stride_in_bytes); + int stbi_write_bmp_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data); + int stbi_write_tga_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data); + int stbi_write_hdr_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const float *data); + int stbi_write_jpg_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data, int quality); + + where the callback is: + void stbi_write_func(void *context, void *data, int size); + + You can configure it with these global variables: + int stbi_write_tga_with_rle; // defaults to true; set to 0 to disable RLE + int stbi_write_png_compression_level; // defaults to 8; set to higher for more compression + int stbi_write_force_png_filter; // defaults to -1; set to 0..5 to force a filter mode + + + You can define STBI_WRITE_NO_STDIO to disable the file variant of these + functions, so the library will not use stdio.h at all. However, this will + also disable HDR writing, because it requires stdio for formatted output. + + Each function returns 0 on failure and non-0 on success. + + The functions create an image file defined by the parameters. The image + is a rectangle of pixels stored from left-to-right, top-to-bottom. + Each pixel contains 'comp' channels of data stored interleaved with 8-bits + per channel, in the following order: 1=Y, 2=YA, 3=RGB, 4=RGBA. (Y is + monochrome color.) The rectangle is 'w' pixels wide and 'h' pixels tall. + The *data pointer points to the first byte of the top-left-most pixel. + For PNG, "stride_in_bytes" is the distance in bytes from the first byte of + a row of pixels to the first byte of the next row of pixels. + + PNG creates output files with the same number of components as the input. + The BMP format expands Y to RGB in the file format and does not + output alpha. + + PNG supports writing rectangles of data even when the bytes storing rows of + data are not consecutive in memory (e.g. sub-rectangles of a larger image), + by supplying the stride between the beginning of adjacent rows. The other + formats do not. (Thus you cannot write a native-format BMP through the BMP + writer, both because it is in BGR order and because it may have padding + at the end of the line.) + + PNG allows you to set the deflate compression level by setting the global + variable 'stbi_write_png_compression_level' (it defaults to 8). + + HDR expects linear float data. Since the format is always 32-bit rgb(e) + data, alpha (if provided) is discarded, and for monochrome data it is + replicated across all three channels. + + TGA supports RLE or non-RLE compressed data. To use non-RLE-compressed + data, set the global variable 'stbi_write_tga_with_rle' to 0. + + JPEG does ignore alpha channels in input data; quality is between 1 and 100. + Higher quality looks better but results in a bigger image. + JPEG baseline (no JPEG progressive). + +CREDITS: + + + Sean Barrett - PNG/BMP/TGA + Baldur Karlsson - HDR + Jean-Sebastien Guay - TGA monochrome + Tim Kelsey - misc enhancements + Alan Hickman - TGA RLE + Emmanuel Julien - initial file IO callback implementation + Jon Olick - original jo_jpeg.cpp code + Daniel Gibson - integrate JPEG, allow external zlib + Aarni Koskela - allow choosing PNG filter + + bugfixes: + github:Chribba + Guillaume Chereau + github:jry2 + github:romigrou + Sergio Gonzalez + Jonas Karlsson + Filip Wasil + Thatcher Ulrich + github:poppolopoppo + Patrick Boettcher + github:xeekworx + Cap Petschulat + Simon Rodriguez + Ivan Tikhonov + github:ignotion + Adam Schackart + +LICENSE + + See end of file for license information. + +*/ + +#ifndef INCLUDE_STB_IMAGE_WRITE_H +#define INCLUDE_STB_IMAGE_WRITE_H + +// if STB_IMAGE_WRITE_STATIC causes problems, try defining STBIWDEF to 'inline' or 'static inline' +#ifndef STBIWDEF +#ifdef STB_IMAGE_WRITE_STATIC +#define STBIWDEF static +#else +#ifdef __cplusplus +#define STBIWDEF extern "C" +#else +#define STBIWDEF extern +#endif +#endif +#endif + +#ifndef STB_IMAGE_WRITE_STATIC // C++ forbids static forward declarations +extern int stbi_write_tga_with_rle; +extern int stbi_write_png_compression_level; +extern int stbi_write_force_png_filter; +#endif + +#ifndef STBI_WRITE_NO_STDIO +STBIWDEF int stbi_write_png(char const *filename, int w, int h, int comp, const void *data, int stride_in_bytes); +STBIWDEF int stbi_write_bmp(char const *filename, int w, int h, int comp, const void *data); +STBIWDEF int stbi_write_tga(char const *filename, int w, int h, int comp, const void *data); +STBIWDEF int stbi_write_hdr(char const *filename, int w, int h, int comp, const float *data); +STBIWDEF int stbi_write_jpg(char const *filename, int x, int y, int comp, const void *data, int quality); +#endif + +typedef void stbi_write_func(void *context, void *data, int size); + +STBIWDEF int stbi_write_png_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data, int stride_in_bytes); +STBIWDEF int stbi_write_bmp_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data); +STBIWDEF int stbi_write_tga_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data); +STBIWDEF int stbi_write_hdr_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const float *data); +STBIWDEF int stbi_write_jpg_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data, int quality); + +STBIWDEF void stbi_flip_vertically_on_write(int flip_boolean); + +#endif//INCLUDE_STB_IMAGE_WRITE_H + +#ifdef STB_IMAGE_WRITE_IMPLEMENTATION + +#ifdef _WIN32 + #ifndef _CRT_SECURE_NO_WARNINGS + #define _CRT_SECURE_NO_WARNINGS + #endif + #ifndef _CRT_NONSTDC_NO_DEPRECATE + #define _CRT_NONSTDC_NO_DEPRECATE + #endif +#endif + +#ifndef STBI_WRITE_NO_STDIO +#include +#endif // STBI_WRITE_NO_STDIO + +#include +#include +#include +#include + +#if defined(STBIW_MALLOC) && defined(STBIW_FREE) && (defined(STBIW_REALLOC) || defined(STBIW_REALLOC_SIZED)) +// ok +#elif !defined(STBIW_MALLOC) && !defined(STBIW_FREE) && !defined(STBIW_REALLOC) && !defined(STBIW_REALLOC_SIZED) +// ok +#else +#error "Must define all or none of STBIW_MALLOC, STBIW_FREE, and STBIW_REALLOC (or STBIW_REALLOC_SIZED)." +#endif + +#ifndef STBIW_MALLOC +#define STBIW_MALLOC(sz) malloc(sz) +#define STBIW_REALLOC(p,newsz) realloc(p,newsz) +#define STBIW_FREE(p) free(p) +#endif + +#ifndef STBIW_REALLOC_SIZED +#define STBIW_REALLOC_SIZED(p,oldsz,newsz) STBIW_REALLOC(p,newsz) +#endif + + +#ifndef STBIW_MEMMOVE +#define STBIW_MEMMOVE(a,b,sz) memmove(a,b,sz) +#endif + + +#ifndef STBIW_ASSERT +#include +#define STBIW_ASSERT(x) assert(x) +#endif + +#define STBIW_UCHAR(x) (unsigned char) ((x) & 0xff) + +#ifdef STB_IMAGE_WRITE_STATIC +static int stbi__flip_vertically_on_write=0; +static int stbi_write_png_compression_level = 8; +static int stbi_write_tga_with_rle = 1; +static int stbi_write_force_png_filter = -1; +#else +int stbi_write_png_compression_level = 8; +int stbi__flip_vertically_on_write=0; +int stbi_write_tga_with_rle = 1; +int stbi_write_force_png_filter = -1; +#endif + +STBIWDEF void stbi_flip_vertically_on_write(int flag) +{ + stbi__flip_vertically_on_write = flag; +} + +typedef struct +{ + stbi_write_func *func; + void *context; +} stbi__write_context; + +// initialize a callback-based context +static void stbi__start_write_callbacks(stbi__write_context *s, stbi_write_func *c, void *context) +{ + s->func = c; + s->context = context; +} + +#ifndef STBI_WRITE_NO_STDIO + +static void stbi__stdio_write(void *context, void *data, int size) +{ + fwrite(data,1,size,(FILE*) context); +} + +static int stbi__start_write_file(stbi__write_context *s, const char *filename) +{ + FILE *f; +#ifdef STBI_MSC_SECURE_CRT + if (fopen_s(&f, filename, "wb")) + f = NULL; +#else + f = fopen(filename, "wb"); +#endif + stbi__start_write_callbacks(s, stbi__stdio_write, (void *) f); + return f != NULL; +} + +static void stbi__end_write_file(stbi__write_context *s) +{ + fclose((FILE *)s->context); +} + +#endif // !STBI_WRITE_NO_STDIO + +typedef unsigned int stbiw_uint32; +typedef int stb_image_write_test[sizeof(stbiw_uint32)==4 ? 1 : -1]; + +static void stbiw__writefv(stbi__write_context *s, const char *fmt, va_list v) +{ + while (*fmt) { + switch (*fmt++) { + case ' ': break; + case '1': { unsigned char x = STBIW_UCHAR(va_arg(v, int)); + s->func(s->context,&x,1); + break; } + case '2': { int x = va_arg(v,int); + unsigned char b[2]; + b[0] = STBIW_UCHAR(x); + b[1] = STBIW_UCHAR(x>>8); + s->func(s->context,b,2); + break; } + case '4': { stbiw_uint32 x = va_arg(v,int); + unsigned char b[4]; + b[0]=STBIW_UCHAR(x); + b[1]=STBIW_UCHAR(x>>8); + b[2]=STBIW_UCHAR(x>>16); + b[3]=STBIW_UCHAR(x>>24); + s->func(s->context,b,4); + break; } + default: + STBIW_ASSERT(0); + return; + } + } +} + +static void stbiw__writef(stbi__write_context *s, const char *fmt, ...) +{ + va_list v; + va_start(v, fmt); + stbiw__writefv(s, fmt, v); + va_end(v); +} + +static void stbiw__putc(stbi__write_context *s, unsigned char c) +{ + s->func(s->context, &c, 1); +} + +static void stbiw__write3(stbi__write_context *s, unsigned char a, unsigned char b, unsigned char c) +{ + unsigned char arr[3]; + arr[0] = a, arr[1] = b, arr[2] = c; + s->func(s->context, arr, 3); +} + +static void stbiw__write_pixel(stbi__write_context *s, int rgb_dir, int comp, int write_alpha, int expand_mono, unsigned char *d) +{ + unsigned char bg[3] = { 255, 0, 255}, px[3]; + int k; + + if (write_alpha < 0) + s->func(s->context, &d[comp - 1], 1); + + switch (comp) { + case 2: // 2 pixels = mono + alpha, alpha is written separately, so same as 1-channel case + case 1: + if (expand_mono) + stbiw__write3(s, d[0], d[0], d[0]); // monochrome bmp + else + s->func(s->context, d, 1); // monochrome TGA + break; + case 4: + if (!write_alpha) { + // composite against pink background + for (k = 0; k < 3; ++k) + px[k] = bg[k] + ((d[k] - bg[k]) * d[3]) / 255; + stbiw__write3(s, px[1 - rgb_dir], px[1], px[1 + rgb_dir]); + break; + } + /* FALLTHROUGH */ + case 3: + stbiw__write3(s, d[1 - rgb_dir], d[1], d[1 + rgb_dir]); + break; + } + if (write_alpha > 0) + s->func(s->context, &d[comp - 1], 1); +} + +static void stbiw__write_pixels(stbi__write_context *s, int rgb_dir, int vdir, int x, int y, int comp, void *data, int write_alpha, int scanline_pad, int expand_mono) +{ + stbiw_uint32 zero = 0; + int i,j, j_end; + + if (y <= 0) + return; + + if (stbi__flip_vertically_on_write) + vdir *= -1; + + if (vdir < 0) + j_end = -1, j = y-1; + else + j_end = y, j = 0; + + for (; j != j_end; j += vdir) { + for (i=0; i < x; ++i) { + unsigned char *d = (unsigned char *) data + (j*x+i)*comp; + stbiw__write_pixel(s, rgb_dir, comp, write_alpha, expand_mono, d); + } + s->func(s->context, &zero, scanline_pad); + } +} + +static int stbiw__outfile(stbi__write_context *s, int rgb_dir, int vdir, int x, int y, int comp, int expand_mono, void *data, int alpha, int pad, const char *fmt, ...) +{ + if (y < 0 || x < 0) { + return 0; + } else { + va_list v; + va_start(v, fmt); + stbiw__writefv(s, fmt, v); + va_end(v); + stbiw__write_pixels(s,rgb_dir,vdir,x,y,comp,data,alpha,pad, expand_mono); + return 1; + } +} + +static int stbi_write_bmp_core(stbi__write_context *s, int x, int y, int comp, const void *data) +{ + int pad = (-x*3) & 3; + return stbiw__outfile(s,-1,-1,x,y,comp,1,(void *) data,0,pad, + "11 4 22 4" "4 44 22 444444", + 'B', 'M', 14+40+(x*3+pad)*y, 0,0, 14+40, // file header + 40, x,y, 1,24, 0,0,0,0,0,0); // bitmap header +} + +STBIWDEF int stbi_write_bmp_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data) +{ + stbi__write_context s; + stbi__start_write_callbacks(&s, func, context); + return stbi_write_bmp_core(&s, x, y, comp, data); +} + +#ifndef STBI_WRITE_NO_STDIO +STBIWDEF int stbi_write_bmp(char const *filename, int x, int y, int comp, const void *data) +{ + stbi__write_context s; + if (stbi__start_write_file(&s,filename)) { + int r = stbi_write_bmp_core(&s, x, y, comp, data); + stbi__end_write_file(&s); + return r; + } else + return 0; +} +#endif //!STBI_WRITE_NO_STDIO + +static int stbi_write_tga_core(stbi__write_context *s, int x, int y, int comp, void *data) +{ + int has_alpha = (comp == 2 || comp == 4); + int colorbytes = has_alpha ? comp-1 : comp; + int format = colorbytes < 2 ? 3 : 2; // 3 color channels (RGB/RGBA) = 2, 1 color channel (Y/YA) = 3 + + if (y < 0 || x < 0) + return 0; + + if (!stbi_write_tga_with_rle) { + return stbiw__outfile(s, -1, -1, x, y, comp, 0, (void *) data, has_alpha, 0, + "111 221 2222 11", 0, 0, format, 0, 0, 0, 0, 0, x, y, (colorbytes + has_alpha) * 8, has_alpha * 8); + } else { + int i,j,k; + int jend, jdir; + + stbiw__writef(s, "111 221 2222 11", 0,0,format+8, 0,0,0, 0,0,x,y, (colorbytes + has_alpha) * 8, has_alpha * 8); + + if (stbi__flip_vertically_on_write) { + j = 0; + jend = y; + jdir = 1; + } else { + j = y-1; + jend = -1; + jdir = -1; + } + for (; j != jend; j += jdir) { + unsigned char *row = (unsigned char *) data + j * x * comp; + int len; + + for (i = 0; i < x; i += len) { + unsigned char *begin = row + i * comp; + int diff = 1; + len = 1; + + if (i < x - 1) { + ++len; + diff = memcmp(begin, row + (i + 1) * comp, comp); + if (diff) { + const unsigned char *prev = begin; + for (k = i + 2; k < x && len < 128; ++k) { + if (memcmp(prev, row + k * comp, comp)) { + prev += comp; + ++len; + } else { + --len; + break; + } + } + } else { + for (k = i + 2; k < x && len < 128; ++k) { + if (!memcmp(begin, row + k * comp, comp)) { + ++len; + } else { + break; + } + } + } + } + + if (diff) { + unsigned char header = STBIW_UCHAR(len - 1); + s->func(s->context, &header, 1); + for (k = 0; k < len; ++k) { + stbiw__write_pixel(s, -1, comp, has_alpha, 0, begin + k * comp); + } + } else { + unsigned char header = STBIW_UCHAR(len - 129); + s->func(s->context, &header, 1); + stbiw__write_pixel(s, -1, comp, has_alpha, 0, begin); + } + } + } + } + return 1; +} + +STBIWDEF int stbi_write_tga_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data) +{ + stbi__write_context s; + stbi__start_write_callbacks(&s, func, context); + return stbi_write_tga_core(&s, x, y, comp, (void *) data); +} + +#ifndef STBI_WRITE_NO_STDIO +STBIWDEF int stbi_write_tga(char const *filename, int x, int y, int comp, const void *data) +{ + stbi__write_context s; + if (stbi__start_write_file(&s,filename)) { + int r = stbi_write_tga_core(&s, x, y, comp, (void *) data); + stbi__end_write_file(&s); + return r; + } else + return 0; +} +#endif + +// ************************************************************************************************* +// Radiance RGBE HDR writer +// by Baldur Karlsson + +#define stbiw__max(a, b) ((a) > (b) ? (a) : (b)) + +void stbiw__linear_to_rgbe(unsigned char *rgbe, float *linear) +{ + int exponent; + float maxcomp = stbiw__max(linear[0], stbiw__max(linear[1], linear[2])); + + if (maxcomp < 1e-32f) { + rgbe[0] = rgbe[1] = rgbe[2] = rgbe[3] = 0; + } else { + float normalize = (float) frexp(maxcomp, &exponent) * 256.0f/maxcomp; + + rgbe[0] = (unsigned char)(linear[0] * normalize); + rgbe[1] = (unsigned char)(linear[1] * normalize); + rgbe[2] = (unsigned char)(linear[2] * normalize); + rgbe[3] = (unsigned char)(exponent + 128); + } +} + +void stbiw__write_run_data(stbi__write_context *s, int length, unsigned char databyte) +{ + unsigned char lengthbyte = STBIW_UCHAR(length+128); + STBIW_ASSERT(length+128 <= 255); + s->func(s->context, &lengthbyte, 1); + s->func(s->context, &databyte, 1); +} + +void stbiw__write_dump_data(stbi__write_context *s, int length, unsigned char *data) +{ + unsigned char lengthbyte = STBIW_UCHAR(length); + STBIW_ASSERT(length <= 128); // inconsistent with spec but consistent with official code + s->func(s->context, &lengthbyte, 1); + s->func(s->context, data, length); +} + +void stbiw__write_hdr_scanline(stbi__write_context *s, int width, int ncomp, unsigned char *scratch, float *scanline) +{ + unsigned char scanlineheader[4] = { 2, 2, 0, 0 }; + unsigned char rgbe[4]; + float linear[3]; + int x; + + scanlineheader[2] = (width&0xff00)>>8; + scanlineheader[3] = (width&0x00ff); + + /* skip RLE for images too small or large */ + if (width < 8 || width >= 32768) { + for (x=0; x < width; x++) { + switch (ncomp) { + case 4: /* fallthrough */ + case 3: linear[2] = scanline[x*ncomp + 2]; + linear[1] = scanline[x*ncomp + 1]; + linear[0] = scanline[x*ncomp + 0]; + break; + default: + linear[0] = linear[1] = linear[2] = scanline[x*ncomp + 0]; + break; + } + stbiw__linear_to_rgbe(rgbe, linear); + s->func(s->context, rgbe, 4); + } + } else { + int c,r; + /* encode into scratch buffer */ + for (x=0; x < width; x++) { + switch(ncomp) { + case 4: /* fallthrough */ + case 3: linear[2] = scanline[x*ncomp + 2]; + linear[1] = scanline[x*ncomp + 1]; + linear[0] = scanline[x*ncomp + 0]; + break; + default: + linear[0] = linear[1] = linear[2] = scanline[x*ncomp + 0]; + break; + } + stbiw__linear_to_rgbe(rgbe, linear); + scratch[x + width*0] = rgbe[0]; + scratch[x + width*1] = rgbe[1]; + scratch[x + width*2] = rgbe[2]; + scratch[x + width*3] = rgbe[3]; + } + + s->func(s->context, scanlineheader, 4); + + /* RLE each component separately */ + for (c=0; c < 4; c++) { + unsigned char *comp = &scratch[width*c]; + + x = 0; + while (x < width) { + // find first run + r = x; + while (r+2 < width) { + if (comp[r] == comp[r+1] && comp[r] == comp[r+2]) + break; + ++r; + } + if (r+2 >= width) + r = width; + // dump up to first run + while (x < r) { + int len = r-x; + if (len > 128) len = 128; + stbiw__write_dump_data(s, len, &comp[x]); + x += len; + } + // if there's a run, output it + if (r+2 < width) { // same test as what we break out of in search loop, so only true if we break'd + // find next byte after run + while (r < width && comp[r] == comp[x]) + ++r; + // output run up to r + while (x < r) { + int len = r-x; + if (len > 127) len = 127; + stbiw__write_run_data(s, len, comp[x]); + x += len; + } + } + } + } + } +} + +static int stbi_write_hdr_core(stbi__write_context *s, int x, int y, int comp, float *data) +{ + if (y <= 0 || x <= 0 || data == NULL) + return 0; + else { + // Each component is stored separately. Allocate scratch space for full output scanline. + unsigned char *scratch = (unsigned char *) STBIW_MALLOC(x*4); + int i, len; + char buffer[128]; + char header[] = "#?RADIANCE\n# Written by stb_image_write.h\nFORMAT=32-bit_rle_rgbe\n"; + s->func(s->context, header, sizeof(header)-1); + +#ifdef STBI_MSC_SECURE_CRT + len = sprintf_s(buffer, "EXPOSURE= 1.0000000000000\n\n-Y %d +X %d\n", y, x); +#else + len = sprintf(buffer, "EXPOSURE= 1.0000000000000\n\n-Y %d +X %d\n", y, x); +#endif + s->func(s->context, buffer, len); + + for(i=0; i < y; i++) + stbiw__write_hdr_scanline(s, x, comp, scratch, data + comp*x*(stbi__flip_vertically_on_write ? y-1-i : i)*x); + STBIW_FREE(scratch); + return 1; + } +} + +STBIWDEF int stbi_write_hdr_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const float *data) +{ + stbi__write_context s; + stbi__start_write_callbacks(&s, func, context); + return stbi_write_hdr_core(&s, x, y, comp, (float *) data); +} + +#ifndef STBI_WRITE_NO_STDIO +STBIWDEF int stbi_write_hdr(char const *filename, int x, int y, int comp, const float *data) +{ + stbi__write_context s; + if (stbi__start_write_file(&s,filename)) { + int r = stbi_write_hdr_core(&s, x, y, comp, (float *) data); + stbi__end_write_file(&s); + return r; + } else + return 0; +} +#endif // STBI_WRITE_NO_STDIO + + +////////////////////////////////////////////////////////////////////////////// +// +// PNG writer +// + +#ifndef STBIW_ZLIB_COMPRESS +// stretchy buffer; stbiw__sbpush() == vector<>::push_back() -- stbiw__sbcount() == vector<>::size() +#define stbiw__sbraw(a) ((int *) (a) - 2) +#define stbiw__sbm(a) stbiw__sbraw(a)[0] +#define stbiw__sbn(a) stbiw__sbraw(a)[1] + +#define stbiw__sbneedgrow(a,n) ((a)==0 || stbiw__sbn(a)+n >= stbiw__sbm(a)) +#define stbiw__sbmaybegrow(a,n) (stbiw__sbneedgrow(a,(n)) ? stbiw__sbgrow(a,n) : 0) +#define stbiw__sbgrow(a,n) stbiw__sbgrowf((void **) &(a), (n), sizeof(*(a))) + +#define stbiw__sbpush(a, v) (stbiw__sbmaybegrow(a,1), (a)[stbiw__sbn(a)++] = (v)) +#define stbiw__sbcount(a) ((a) ? stbiw__sbn(a) : 0) +#define stbiw__sbfree(a) ((a) ? STBIW_FREE(stbiw__sbraw(a)),0 : 0) + +static void *stbiw__sbgrowf(void **arr, int increment, int itemsize) +{ + int m = *arr ? 2*stbiw__sbm(*arr)+increment : increment+1; + void *p = STBIW_REALLOC_SIZED(*arr ? stbiw__sbraw(*arr) : 0, *arr ? (stbiw__sbm(*arr)*itemsize + sizeof(int)*2) : 0, itemsize * m + sizeof(int)*2); + STBIW_ASSERT(p); + if (p) { + if (!*arr) ((int *) p)[1] = 0; + *arr = (void *) ((int *) p + 2); + stbiw__sbm(*arr) = m; + } + return *arr; +} + +static unsigned char *stbiw__zlib_flushf(unsigned char *data, unsigned int *bitbuffer, int *bitcount) +{ + while (*bitcount >= 8) { + stbiw__sbpush(data, STBIW_UCHAR(*bitbuffer)); + *bitbuffer >>= 8; + *bitcount -= 8; + } + return data; +} + +static int stbiw__zlib_bitrev(int code, int codebits) +{ + int res=0; + while (codebits--) { + res = (res << 1) | (code & 1); + code >>= 1; + } + return res; +} + +static unsigned int stbiw__zlib_countm(unsigned char *a, unsigned char *b, int limit) +{ + int i; + for (i=0; i < limit && i < 258; ++i) + if (a[i] != b[i]) break; + return i; +} + +static unsigned int stbiw__zhash(unsigned char *data) +{ + stbiw_uint32 hash = data[0] + (data[1] << 8) + (data[2] << 16); + hash ^= hash << 3; + hash += hash >> 5; + hash ^= hash << 4; + hash += hash >> 17; + hash ^= hash << 25; + hash += hash >> 6; + return hash; +} + +#define stbiw__zlib_flush() (out = stbiw__zlib_flushf(out, &bitbuf, &bitcount)) +#define stbiw__zlib_add(code,codebits) \ + (bitbuf |= (code) << bitcount, bitcount += (codebits), stbiw__zlib_flush()) +#define stbiw__zlib_huffa(b,c) stbiw__zlib_add(stbiw__zlib_bitrev(b,c),c) +// default huffman tables +#define stbiw__zlib_huff1(n) stbiw__zlib_huffa(0x30 + (n), 8) +#define stbiw__zlib_huff2(n) stbiw__zlib_huffa(0x190 + (n)-144, 9) +#define stbiw__zlib_huff3(n) stbiw__zlib_huffa(0 + (n)-256,7) +#define stbiw__zlib_huff4(n) stbiw__zlib_huffa(0xc0 + (n)-280,8) +#define stbiw__zlib_huff(n) ((n) <= 143 ? stbiw__zlib_huff1(n) : (n) <= 255 ? stbiw__zlib_huff2(n) : (n) <= 279 ? stbiw__zlib_huff3(n) : stbiw__zlib_huff4(n)) +#define stbiw__zlib_huffb(n) ((n) <= 143 ? stbiw__zlib_huff1(n) : stbiw__zlib_huff2(n)) + +#define stbiw__ZHASH 16384 + +#endif // STBIW_ZLIB_COMPRESS + +unsigned char * stbi_zlib_compress(unsigned char *data, int data_len, int *out_len, int quality) +{ +#ifdef STBIW_ZLIB_COMPRESS + // user provided a zlib compress implementation, use that + return STBIW_ZLIB_COMPRESS(data, data_len, out_len, quality); +#else // use builtin + static unsigned short lengthc[] = { 3,4,5,6,7,8,9,10,11,13,15,17,19,23,27,31,35,43,51,59,67,83,99,115,131,163,195,227,258, 259 }; + static unsigned char lengtheb[]= { 0,0,0,0,0,0,0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 0 }; + static unsigned short distc[] = { 1,2,3,4,5,7,9,13,17,25,33,49,65,97,129,193,257,385,513,769,1025,1537,2049,3073,4097,6145,8193,12289,16385,24577, 32768 }; + static unsigned char disteb[] = { 0,0,0,0,1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,11,11,12,12,13,13 }; + unsigned int bitbuf=0; + int i,j, bitcount=0; + unsigned char *out = NULL; + unsigned char ***hash_table = (unsigned char***) STBIW_MALLOC(stbiw__ZHASH * sizeof(char**)); + if (hash_table == NULL) + return NULL; + if (quality < 5) quality = 5; + + stbiw__sbpush(out, 0x78); // DEFLATE 32K window + stbiw__sbpush(out, 0x5e); // FLEVEL = 1 + stbiw__zlib_add(1,1); // BFINAL = 1 + stbiw__zlib_add(1,2); // BTYPE = 1 -- fixed huffman + + for (i=0; i < stbiw__ZHASH; ++i) + hash_table[i] = NULL; + + i=0; + while (i < data_len-3) { + // hash next 3 bytes of data to be compressed + int h = stbiw__zhash(data+i)&(stbiw__ZHASH-1), best=3; + unsigned char *bestloc = 0; + unsigned char **hlist = hash_table[h]; + int n = stbiw__sbcount(hlist); + for (j=0; j < n; ++j) { + if (hlist[j]-data > i-32768) { // if entry lies within window + int d = stbiw__zlib_countm(hlist[j], data+i, data_len-i); + if (d >= best) best=d,bestloc=hlist[j]; + } + } + // when hash table entry is too long, delete half the entries + if (hash_table[h] && stbiw__sbn(hash_table[h]) == 2*quality) { + STBIW_MEMMOVE(hash_table[h], hash_table[h]+quality, sizeof(hash_table[h][0])*quality); + stbiw__sbn(hash_table[h]) = quality; + } + stbiw__sbpush(hash_table[h],data+i); + + if (bestloc) { + // "lazy matching" - check match at *next* byte, and if it's better, do cur byte as literal + h = stbiw__zhash(data+i+1)&(stbiw__ZHASH-1); + hlist = hash_table[h]; + n = stbiw__sbcount(hlist); + for (j=0; j < n; ++j) { + if (hlist[j]-data > i-32767) { + int e = stbiw__zlib_countm(hlist[j], data+i+1, data_len-i-1); + if (e > best) { // if next match is better, bail on current match + bestloc = NULL; + break; + } + } + } + } + + if (bestloc) { + int d = (int) (data+i - bestloc); // distance back + STBIW_ASSERT(d <= 32767 && best <= 258); + for (j=0; best > lengthc[j+1]-1; ++j); + stbiw__zlib_huff(j+257); + if (lengtheb[j]) stbiw__zlib_add(best - lengthc[j], lengtheb[j]); + for (j=0; d > distc[j+1]-1; ++j); + stbiw__zlib_add(stbiw__zlib_bitrev(j,5),5); + if (disteb[j]) stbiw__zlib_add(d - distc[j], disteb[j]); + i += best; + } else { + stbiw__zlib_huffb(data[i]); + ++i; + } + } + // write out final bytes + for (;i < data_len; ++i) + stbiw__zlib_huffb(data[i]); + stbiw__zlib_huff(256); // end of block + // pad with 0 bits to byte boundary + while (bitcount) + stbiw__zlib_add(0,1); + + for (i=0; i < stbiw__ZHASH; ++i) + (void) stbiw__sbfree(hash_table[i]); + STBIW_FREE(hash_table); + + { + // compute adler32 on input + unsigned int s1=1, s2=0; + int blocklen = (int) (data_len % 5552); + j=0; + while (j < data_len) { + for (i=0; i < blocklen; ++i) s1 += data[j+i], s2 += s1; + s1 %= 65521, s2 %= 65521; + j += blocklen; + blocklen = 5552; + } + stbiw__sbpush(out, STBIW_UCHAR(s2 >> 8)); + stbiw__sbpush(out, STBIW_UCHAR(s2)); + stbiw__sbpush(out, STBIW_UCHAR(s1 >> 8)); + stbiw__sbpush(out, STBIW_UCHAR(s1)); + } + *out_len = stbiw__sbn(out); + // make returned pointer freeable + STBIW_MEMMOVE(stbiw__sbraw(out), out, *out_len); + return (unsigned char *) stbiw__sbraw(out); +#endif // STBIW_ZLIB_COMPRESS +} + +static unsigned int stbiw__crc32(unsigned char *buffer, int len) +{ + static unsigned int crc_table[256] = + { + 0x00000000, 0x77073096, 0xEE0E612C, 0x990951BA, 0x076DC419, 0x706AF48F, 0xE963A535, 0x9E6495A3, + 0x0eDB8832, 0x79DCB8A4, 0xE0D5E91E, 0x97D2D988, 0x09B64C2B, 0x7EB17CBD, 0xE7B82D07, 0x90BF1D91, + 0x1DB71064, 0x6AB020F2, 0xF3B97148, 0x84BE41DE, 0x1ADAD47D, 0x6DDDE4EB, 0xF4D4B551, 0x83D385C7, + 0x136C9856, 0x646BA8C0, 0xFD62F97A, 0x8A65C9EC, 0x14015C4F, 0x63066CD9, 0xFA0F3D63, 0x8D080DF5, + 0x3B6E20C8, 0x4C69105E, 0xD56041E4, 0xA2677172, 0x3C03E4D1, 0x4B04D447, 0xD20D85FD, 0xA50AB56B, + 0x35B5A8FA, 0x42B2986C, 0xDBBBC9D6, 0xACBCF940, 0x32D86CE3, 0x45DF5C75, 0xDCD60DCF, 0xABD13D59, + 0x26D930AC, 0x51DE003A, 0xC8D75180, 0xBFD06116, 0x21B4F4B5, 0x56B3C423, 0xCFBA9599, 0xB8BDA50F, + 0x2802B89E, 0x5F058808, 0xC60CD9B2, 0xB10BE924, 0x2F6F7C87, 0x58684C11, 0xC1611DAB, 0xB6662D3D, + 0x76DC4190, 0x01DB7106, 0x98D220BC, 0xEFD5102A, 0x71B18589, 0x06B6B51F, 0x9FBFE4A5, 0xE8B8D433, + 0x7807C9A2, 0x0F00F934, 0x9609A88E, 0xE10E9818, 0x7F6A0DBB, 0x086D3D2D, 0x91646C97, 0xE6635C01, + 0x6B6B51F4, 0x1C6C6162, 0x856530D8, 0xF262004E, 0x6C0695ED, 0x1B01A57B, 0x8208F4C1, 0xF50FC457, + 0x65B0D9C6, 0x12B7E950, 0x8BBEB8EA, 0xFCB9887C, 0x62DD1DDF, 0x15DA2D49, 0x8CD37CF3, 0xFBD44C65, + 0x4DB26158, 0x3AB551CE, 0xA3BC0074, 0xD4BB30E2, 0x4ADFA541, 0x3DD895D7, 0xA4D1C46D, 0xD3D6F4FB, + 0x4369E96A, 0x346ED9FC, 0xAD678846, 0xDA60B8D0, 0x44042D73, 0x33031DE5, 0xAA0A4C5F, 0xDD0D7CC9, + 0x5005713C, 0x270241AA, 0xBE0B1010, 0xC90C2086, 0x5768B525, 0x206F85B3, 0xB966D409, 0xCE61E49F, + 0x5EDEF90E, 0x29D9C998, 0xB0D09822, 0xC7D7A8B4, 0x59B33D17, 0x2EB40D81, 0xB7BD5C3B, 0xC0BA6CAD, + 0xEDB88320, 0x9ABFB3B6, 0x03B6E20C, 0x74B1D29A, 0xEAD54739, 0x9DD277AF, 0x04DB2615, 0x73DC1683, + 0xE3630B12, 0x94643B84, 0x0D6D6A3E, 0x7A6A5AA8, 0xE40ECF0B, 0x9309FF9D, 0x0A00AE27, 0x7D079EB1, + 0xF00F9344, 0x8708A3D2, 0x1E01F268, 0x6906C2FE, 0xF762575D, 0x806567CB, 0x196C3671, 0x6E6B06E7, + 0xFED41B76, 0x89D32BE0, 0x10DA7A5A, 0x67DD4ACC, 0xF9B9DF6F, 0x8EBEEFF9, 0x17B7BE43, 0x60B08ED5, + 0xD6D6A3E8, 0xA1D1937E, 0x38D8C2C4, 0x4FDFF252, 0xD1BB67F1, 0xA6BC5767, 0x3FB506DD, 0x48B2364B, + 0xD80D2BDA, 0xAF0A1B4C, 0x36034AF6, 0x41047A60, 0xDF60EFC3, 0xA867DF55, 0x316E8EEF, 0x4669BE79, + 0xCB61B38C, 0xBC66831A, 0x256FD2A0, 0x5268E236, 0xCC0C7795, 0xBB0B4703, 0x220216B9, 0x5505262F, + 0xC5BA3BBE, 0xB2BD0B28, 0x2BB45A92, 0x5CB36A04, 0xC2D7FFA7, 0xB5D0CF31, 0x2CD99E8B, 0x5BDEAE1D, + 0x9B64C2B0, 0xEC63F226, 0x756AA39C, 0x026D930A, 0x9C0906A9, 0xEB0E363F, 0x72076785, 0x05005713, + 0x95BF4A82, 0xE2B87A14, 0x7BB12BAE, 0x0CB61B38, 0x92D28E9B, 0xE5D5BE0D, 0x7CDCEFB7, 0x0BDBDF21, + 0x86D3D2D4, 0xF1D4E242, 0x68DDB3F8, 0x1FDA836E, 0x81BE16CD, 0xF6B9265B, 0x6FB077E1, 0x18B74777, + 0x88085AE6, 0xFF0F6A70, 0x66063BCA, 0x11010B5C, 0x8F659EFF, 0xF862AE69, 0x616BFFD3, 0x166CCF45, + 0xA00AE278, 0xD70DD2EE, 0x4E048354, 0x3903B3C2, 0xA7672661, 0xD06016F7, 0x4969474D, 0x3E6E77DB, + 0xAED16A4A, 0xD9D65ADC, 0x40DF0B66, 0x37D83BF0, 0xA9BCAE53, 0xDEBB9EC5, 0x47B2CF7F, 0x30B5FFE9, + 0xBDBDF21C, 0xCABAC28A, 0x53B39330, 0x24B4A3A6, 0xBAD03605, 0xCDD70693, 0x54DE5729, 0x23D967BF, + 0xB3667A2E, 0xC4614AB8, 0x5D681B02, 0x2A6F2B94, 0xB40BBE37, 0xC30C8EA1, 0x5A05DF1B, 0x2D02EF8D + }; + + unsigned int crc = ~0u; + int i; + for (i=0; i < len; ++i) + crc = (crc >> 8) ^ crc_table[buffer[i] ^ (crc & 0xff)]; + return ~crc; +} + +#define stbiw__wpng4(o,a,b,c,d) ((o)[0]=STBIW_UCHAR(a),(o)[1]=STBIW_UCHAR(b),(o)[2]=STBIW_UCHAR(c),(o)[3]=STBIW_UCHAR(d),(o)+=4) +#define stbiw__wp32(data,v) stbiw__wpng4(data, (v)>>24,(v)>>16,(v)>>8,(v)); +#define stbiw__wptag(data,s) stbiw__wpng4(data, s[0],s[1],s[2],s[3]) + +static void stbiw__wpcrc(unsigned char **data, int len) +{ + unsigned int crc = stbiw__crc32(*data - len - 4, len+4); + stbiw__wp32(*data, crc); +} + +static unsigned char stbiw__paeth(int a, int b, int c) +{ + int p = a + b - c, pa = abs(p-a), pb = abs(p-b), pc = abs(p-c); + if (pa <= pb && pa <= pc) return STBIW_UCHAR(a); + if (pb <= pc) return STBIW_UCHAR(b); + return STBIW_UCHAR(c); +} + +// @OPTIMIZE: provide an option that always forces left-predict or paeth predict +static void stbiw__encode_png_line(unsigned char *pixels, int stride_bytes, int width, int height, int y, int n, int filter_type, signed char *line_buffer) +{ + static int mapping[] = { 0,1,2,3,4 }; + static int firstmap[] = { 0,1,0,5,6 }; + int *mymap = (y != 0) ? mapping : firstmap; + int i; + int type = mymap[filter_type]; + unsigned char *z = pixels + stride_bytes * (stbi__flip_vertically_on_write ? height-1-y : y); + int signed_stride = stbi__flip_vertically_on_write ? -stride_bytes : stride_bytes; + for (i = 0; i < n; ++i) { + switch (type) { + case 0: line_buffer[i] = z[i]; break; + case 1: line_buffer[i] = z[i]; break; + case 2: line_buffer[i] = z[i] - z[i-signed_stride]; break; + case 3: line_buffer[i] = z[i] - (z[i-signed_stride]>>1); break; + case 4: line_buffer[i] = (signed char) (z[i] - stbiw__paeth(0,z[i-signed_stride],0)); break; + case 5: line_buffer[i] = z[i]; break; + case 6: line_buffer[i] = z[i]; break; + } + } + for (i=n; i < width*n; ++i) { + switch (type) { + case 0: line_buffer[i] = z[i]; break; + case 1: line_buffer[i] = z[i] - z[i-n]; break; + case 2: line_buffer[i] = z[i] - z[i-signed_stride]; break; + case 3: line_buffer[i] = z[i] - ((z[i-n] + z[i-signed_stride])>>1); break; + case 4: line_buffer[i] = z[i] - stbiw__paeth(z[i-n], z[i-signed_stride], z[i-signed_stride-n]); break; + case 5: line_buffer[i] = z[i] - (z[i-n]>>1); break; + case 6: line_buffer[i] = z[i] - stbiw__paeth(z[i-n], 0,0); break; + } + } +} + +unsigned char *stbi_write_png_to_mem(unsigned char *pixels, int stride_bytes, int x, int y, int n, int *out_len) +{ + int force_filter = stbi_write_force_png_filter; + int ctype[5] = { -1, 0, 4, 2, 6 }; + unsigned char sig[8] = { 137,80,78,71,13,10,26,10 }; + unsigned char *out,*o, *filt, *zlib; + signed char *line_buffer; + int j,zlen; + + if (stride_bytes == 0) + stride_bytes = x * n; + + if (force_filter >= 5) { + force_filter = -1; + } + + filt = (unsigned char *) STBIW_MALLOC((x*n+1) * y); if (!filt) return 0; + line_buffer = (signed char *) STBIW_MALLOC(x * n); if (!line_buffer) { STBIW_FREE(filt); return 0; } + for (j=0; j < y; ++j) { + int filter_type; + if (force_filter > -1) { + filter_type = force_filter; + stbiw__encode_png_line(pixels, stride_bytes, x, y, j, n, force_filter, line_buffer); + } else { // Estimate the best filter by running through all of them: + int best_filter = 0, best_filter_val = 0x7fffffff, est, i; + for (filter_type = 0; filter_type < 5; filter_type++) { + stbiw__encode_png_line(pixels, stride_bytes, x, y, j, n, filter_type, line_buffer); + + // Estimate the entropy of the line using this filter; the less, the better. + est = 0; + for (i = 0; i < x*n; ++i) { + est += abs((signed char) line_buffer[i]); + } + if (est < best_filter_val) { + best_filter_val = est; + best_filter = filter_type; + } + } + if (filter_type != best_filter) { // If the last iteration already got us the best filter, don't redo it + stbiw__encode_png_line(pixels, stride_bytes, x, y, j, n, best_filter, line_buffer); + filter_type = best_filter; + } + } + // when we get here, filter_type contains the filter type, and line_buffer contains the data + filt[j*(x*n+1)] = (unsigned char) filter_type; + STBIW_MEMMOVE(filt+j*(x*n+1)+1, line_buffer, x*n); + } + STBIW_FREE(line_buffer); + zlib = stbi_zlib_compress(filt, y*( x*n+1), &zlen, stbi_write_png_compression_level); + STBIW_FREE(filt); + if (!zlib) return 0; + + // each tag requires 12 bytes of overhead + out = (unsigned char *) STBIW_MALLOC(8 + 12+13 + 12+zlen + 12); + if (!out) return 0; + *out_len = 8 + 12+13 + 12+zlen + 12; + + o=out; + STBIW_MEMMOVE(o,sig,8); o+= 8; + stbiw__wp32(o, 13); // header length + stbiw__wptag(o, "IHDR"); + stbiw__wp32(o, x); + stbiw__wp32(o, y); + *o++ = 8; + *o++ = STBIW_UCHAR(ctype[n]); + *o++ = 0; + *o++ = 0; + *o++ = 0; + stbiw__wpcrc(&o,13); + + stbiw__wp32(o, zlen); + stbiw__wptag(o, "IDAT"); + STBIW_MEMMOVE(o, zlib, zlen); + o += zlen; + STBIW_FREE(zlib); + stbiw__wpcrc(&o, zlen); + + stbiw__wp32(o,0); + stbiw__wptag(o, "IEND"); + stbiw__wpcrc(&o,0); + + STBIW_ASSERT(o == out + *out_len); + + return out; +} + +#ifndef STBI_WRITE_NO_STDIO +STBIWDEF int stbi_write_png(char const *filename, int x, int y, int comp, const void *data, int stride_bytes) +{ + FILE *f; + int len; + unsigned char *png = stbi_write_png_to_mem((unsigned char *) data, stride_bytes, x, y, comp, &len); + if (png == NULL) return 0; +#ifdef STBI_MSC_SECURE_CRT + if (fopen_s(&f, filename, "wb")) + f = NULL; +#else + f = fopen(filename, "wb"); +#endif + if (!f) { STBIW_FREE(png); return 0; } + fwrite(png, 1, len, f); + fclose(f); + STBIW_FREE(png); + return 1; +} +#endif + +STBIWDEF int stbi_write_png_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data, int stride_bytes) +{ + int len; + unsigned char *png = stbi_write_png_to_mem((unsigned char *) data, stride_bytes, x, y, comp, &len); + if (png == NULL) return 0; + func(context, png, len); + STBIW_FREE(png); + return 1; +} + + +/* *************************************************************************** + * + * JPEG writer + * + * This is based on Jon Olick's jo_jpeg.cpp: + * public domain Simple, Minimalistic JPEG writer - http://www.jonolick.com/code.html + */ + +static const unsigned char stbiw__jpg_ZigZag[] = { 0,1,5,6,14,15,27,28,2,4,7,13,16,26,29,42,3,8,12,17,25,30,41,43,9,11,18, + 24,31,40,44,53,10,19,23,32,39,45,52,54,20,22,33,38,46,51,55,60,21,34,37,47,50,56,59,61,35,36,48,49,57,58,62,63 }; + +static void stbiw__jpg_writeBits(stbi__write_context *s, int *bitBufP, int *bitCntP, const unsigned short *bs) { + int bitBuf = *bitBufP, bitCnt = *bitCntP; + bitCnt += bs[1]; + bitBuf |= bs[0] << (24 - bitCnt); + while(bitCnt >= 8) { + unsigned char c = (bitBuf >> 16) & 255; + stbiw__putc(s, c); + if(c == 255) { + stbiw__putc(s, 0); + } + bitBuf <<= 8; + bitCnt -= 8; + } + *bitBufP = bitBuf; + *bitCntP = bitCnt; +} + +static void stbiw__jpg_DCT(float *d0p, float *d1p, float *d2p, float *d3p, float *d4p, float *d5p, float *d6p, float *d7p) { + float d0 = *d0p, d1 = *d1p, d2 = *d2p, d3 = *d3p, d4 = *d4p, d5 = *d5p, d6 = *d6p, d7 = *d7p; + float z1, z2, z3, z4, z5, z11, z13; + + float tmp0 = d0 + d7; + float tmp7 = d0 - d7; + float tmp1 = d1 + d6; + float tmp6 = d1 - d6; + float tmp2 = d2 + d5; + float tmp5 = d2 - d5; + float tmp3 = d3 + d4; + float tmp4 = d3 - d4; + + // Even part + float tmp10 = tmp0 + tmp3; // phase 2 + float tmp13 = tmp0 - tmp3; + float tmp11 = tmp1 + tmp2; + float tmp12 = tmp1 - tmp2; + + d0 = tmp10 + tmp11; // phase 3 + d4 = tmp10 - tmp11; + + z1 = (tmp12 + tmp13) * 0.707106781f; // c4 + d2 = tmp13 + z1; // phase 5 + d6 = tmp13 - z1; + + // Odd part + tmp10 = tmp4 + tmp5; // phase 2 + tmp11 = tmp5 + tmp6; + tmp12 = tmp6 + tmp7; + + // The rotator is modified from fig 4-8 to avoid extra negations. + z5 = (tmp10 - tmp12) * 0.382683433f; // c6 + z2 = tmp10 * 0.541196100f + z5; // c2-c6 + z4 = tmp12 * 1.306562965f + z5; // c2+c6 + z3 = tmp11 * 0.707106781f; // c4 + + z11 = tmp7 + z3; // phase 5 + z13 = tmp7 - z3; + + *d5p = z13 + z2; // phase 6 + *d3p = z13 - z2; + *d1p = z11 + z4; + *d7p = z11 - z4; + + *d0p = d0; *d2p = d2; *d4p = d4; *d6p = d6; +} + +static void stbiw__jpg_calcBits(int val, unsigned short bits[2]) { + int tmp1 = val < 0 ? -val : val; + val = val < 0 ? val-1 : val; + bits[1] = 1; + while(tmp1 >>= 1) { + ++bits[1]; + } + bits[0] = val & ((1<0)&&(DU[end0pos]==0); --end0pos) { + } + // end0pos = first element in reverse order !=0 + if(end0pos == 0) { + stbiw__jpg_writeBits(s, bitBuf, bitCnt, EOB); + return DU[0]; + } + for(i = 1; i <= end0pos; ++i) { + int startpos = i; + int nrzeroes; + unsigned short bits[2]; + for (; DU[i]==0 && i<=end0pos; ++i) { + } + nrzeroes = i-startpos; + if ( nrzeroes >= 16 ) { + int lng = nrzeroes>>4; + int nrmarker; + for (nrmarker=1; nrmarker <= lng; ++nrmarker) + stbiw__jpg_writeBits(s, bitBuf, bitCnt, M16zeroes); + nrzeroes &= 15; + } + stbiw__jpg_calcBits(DU[i], bits); + stbiw__jpg_writeBits(s, bitBuf, bitCnt, HTAC[(nrzeroes<<4)+bits[1]]); + stbiw__jpg_writeBits(s, bitBuf, bitCnt, bits); + } + if(end0pos != 63) { + stbiw__jpg_writeBits(s, bitBuf, bitCnt, EOB); + } + return DU[0]; +} + +static int stbi_write_jpg_core(stbi__write_context *s, int width, int height, int comp, const void* data, int quality) { + // Constants that don't pollute global namespace + static const unsigned char std_dc_luminance_nrcodes[] = {0,0,1,5,1,1,1,1,1,1,0,0,0,0,0,0,0}; + static const unsigned char std_dc_luminance_values[] = {0,1,2,3,4,5,6,7,8,9,10,11}; + static const unsigned char std_ac_luminance_nrcodes[] = {0,0,2,1,3,3,2,4,3,5,5,4,4,0,0,1,0x7d}; + static const unsigned char std_ac_luminance_values[] = { + 0x01,0x02,0x03,0x00,0x04,0x11,0x05,0x12,0x21,0x31,0x41,0x06,0x13,0x51,0x61,0x07,0x22,0x71,0x14,0x32,0x81,0x91,0xa1,0x08, + 0x23,0x42,0xb1,0xc1,0x15,0x52,0xd1,0xf0,0x24,0x33,0x62,0x72,0x82,0x09,0x0a,0x16,0x17,0x18,0x19,0x1a,0x25,0x26,0x27,0x28, + 0x29,0x2a,0x34,0x35,0x36,0x37,0x38,0x39,0x3a,0x43,0x44,0x45,0x46,0x47,0x48,0x49,0x4a,0x53,0x54,0x55,0x56,0x57,0x58,0x59, + 0x5a,0x63,0x64,0x65,0x66,0x67,0x68,0x69,0x6a,0x73,0x74,0x75,0x76,0x77,0x78,0x79,0x7a,0x83,0x84,0x85,0x86,0x87,0x88,0x89, + 0x8a,0x92,0x93,0x94,0x95,0x96,0x97,0x98,0x99,0x9a,0xa2,0xa3,0xa4,0xa5,0xa6,0xa7,0xa8,0xa9,0xaa,0xb2,0xb3,0xb4,0xb5,0xb6, + 0xb7,0xb8,0xb9,0xba,0xc2,0xc3,0xc4,0xc5,0xc6,0xc7,0xc8,0xc9,0xca,0xd2,0xd3,0xd4,0xd5,0xd6,0xd7,0xd8,0xd9,0xda,0xe1,0xe2, + 0xe3,0xe4,0xe5,0xe6,0xe7,0xe8,0xe9,0xea,0xf1,0xf2,0xf3,0xf4,0xf5,0xf6,0xf7,0xf8,0xf9,0xfa + }; + static const unsigned char std_dc_chrominance_nrcodes[] = {0,0,3,1,1,1,1,1,1,1,1,1,0,0,0,0,0}; + static const unsigned char std_dc_chrominance_values[] = {0,1,2,3,4,5,6,7,8,9,10,11}; + static const unsigned char std_ac_chrominance_nrcodes[] = {0,0,2,1,2,4,4,3,4,7,5,4,4,0,1,2,0x77}; + static const unsigned char std_ac_chrominance_values[] = { + 0x00,0x01,0x02,0x03,0x11,0x04,0x05,0x21,0x31,0x06,0x12,0x41,0x51,0x07,0x61,0x71,0x13,0x22,0x32,0x81,0x08,0x14,0x42,0x91, + 0xa1,0xb1,0xc1,0x09,0x23,0x33,0x52,0xf0,0x15,0x62,0x72,0xd1,0x0a,0x16,0x24,0x34,0xe1,0x25,0xf1,0x17,0x18,0x19,0x1a,0x26, + 0x27,0x28,0x29,0x2a,0x35,0x36,0x37,0x38,0x39,0x3a,0x43,0x44,0x45,0x46,0x47,0x48,0x49,0x4a,0x53,0x54,0x55,0x56,0x57,0x58, + 0x59,0x5a,0x63,0x64,0x65,0x66,0x67,0x68,0x69,0x6a,0x73,0x74,0x75,0x76,0x77,0x78,0x79,0x7a,0x82,0x83,0x84,0x85,0x86,0x87, + 0x88,0x89,0x8a,0x92,0x93,0x94,0x95,0x96,0x97,0x98,0x99,0x9a,0xa2,0xa3,0xa4,0xa5,0xa6,0xa7,0xa8,0xa9,0xaa,0xb2,0xb3,0xb4, + 0xb5,0xb6,0xb7,0xb8,0xb9,0xba,0xc2,0xc3,0xc4,0xc5,0xc6,0xc7,0xc8,0xc9,0xca,0xd2,0xd3,0xd4,0xd5,0xd6,0xd7,0xd8,0xd9,0xda, + 0xe2,0xe3,0xe4,0xe5,0xe6,0xe7,0xe8,0xe9,0xea,0xf2,0xf3,0xf4,0xf5,0xf6,0xf7,0xf8,0xf9,0xfa + }; + // Huffman tables + static const unsigned short YDC_HT[256][2] = { {0,2},{2,3},{3,3},{4,3},{5,3},{6,3},{14,4},{30,5},{62,6},{126,7},{254,8},{510,9}}; + static const unsigned short UVDC_HT[256][2] = { {0,2},{1,2},{2,2},{6,3},{14,4},{30,5},{62,6},{126,7},{254,8},{510,9},{1022,10},{2046,11}}; + static const unsigned short YAC_HT[256][2] = { + {10,4},{0,2},{1,2},{4,3},{11,4},{26,5},{120,7},{248,8},{1014,10},{65410,16},{65411,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {12,4},{27,5},{121,7},{502,9},{2038,11},{65412,16},{65413,16},{65414,16},{65415,16},{65416,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {28,5},{249,8},{1015,10},{4084,12},{65417,16},{65418,16},{65419,16},{65420,16},{65421,16},{65422,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {58,6},{503,9},{4085,12},{65423,16},{65424,16},{65425,16},{65426,16},{65427,16},{65428,16},{65429,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {59,6},{1016,10},{65430,16},{65431,16},{65432,16},{65433,16},{65434,16},{65435,16},{65436,16},{65437,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {122,7},{2039,11},{65438,16},{65439,16},{65440,16},{65441,16},{65442,16},{65443,16},{65444,16},{65445,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {123,7},{4086,12},{65446,16},{65447,16},{65448,16},{65449,16},{65450,16},{65451,16},{65452,16},{65453,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {250,8},{4087,12},{65454,16},{65455,16},{65456,16},{65457,16},{65458,16},{65459,16},{65460,16},{65461,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {504,9},{32704,15},{65462,16},{65463,16},{65464,16},{65465,16},{65466,16},{65467,16},{65468,16},{65469,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {505,9},{65470,16},{65471,16},{65472,16},{65473,16},{65474,16},{65475,16},{65476,16},{65477,16},{65478,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {506,9},{65479,16},{65480,16},{65481,16},{65482,16},{65483,16},{65484,16},{65485,16},{65486,16},{65487,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {1017,10},{65488,16},{65489,16},{65490,16},{65491,16},{65492,16},{65493,16},{65494,16},{65495,16},{65496,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {1018,10},{65497,16},{65498,16},{65499,16},{65500,16},{65501,16},{65502,16},{65503,16},{65504,16},{65505,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {2040,11},{65506,16},{65507,16},{65508,16},{65509,16},{65510,16},{65511,16},{65512,16},{65513,16},{65514,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {65515,16},{65516,16},{65517,16},{65518,16},{65519,16},{65520,16},{65521,16},{65522,16},{65523,16},{65524,16},{0,0},{0,0},{0,0},{0,0},{0,0}, + {2041,11},{65525,16},{65526,16},{65527,16},{65528,16},{65529,16},{65530,16},{65531,16},{65532,16},{65533,16},{65534,16},{0,0},{0,0},{0,0},{0,0},{0,0} + }; + static const unsigned short UVAC_HT[256][2] = { + {0,2},{1,2},{4,3},{10,4},{24,5},{25,5},{56,6},{120,7},{500,9},{1014,10},{4084,12},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {11,4},{57,6},{246,8},{501,9},{2038,11},{4085,12},{65416,16},{65417,16},{65418,16},{65419,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {26,5},{247,8},{1015,10},{4086,12},{32706,15},{65420,16},{65421,16},{65422,16},{65423,16},{65424,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {27,5},{248,8},{1016,10},{4087,12},{65425,16},{65426,16},{65427,16},{65428,16},{65429,16},{65430,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {58,6},{502,9},{65431,16},{65432,16},{65433,16},{65434,16},{65435,16},{65436,16},{65437,16},{65438,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {59,6},{1017,10},{65439,16},{65440,16},{65441,16},{65442,16},{65443,16},{65444,16},{65445,16},{65446,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {121,7},{2039,11},{65447,16},{65448,16},{65449,16},{65450,16},{65451,16},{65452,16},{65453,16},{65454,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {122,7},{2040,11},{65455,16},{65456,16},{65457,16},{65458,16},{65459,16},{65460,16},{65461,16},{65462,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {249,8},{65463,16},{65464,16},{65465,16},{65466,16},{65467,16},{65468,16},{65469,16},{65470,16},{65471,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {503,9},{65472,16},{65473,16},{65474,16},{65475,16},{65476,16},{65477,16},{65478,16},{65479,16},{65480,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {504,9},{65481,16},{65482,16},{65483,16},{65484,16},{65485,16},{65486,16},{65487,16},{65488,16},{65489,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {505,9},{65490,16},{65491,16},{65492,16},{65493,16},{65494,16},{65495,16},{65496,16},{65497,16},{65498,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {506,9},{65499,16},{65500,16},{65501,16},{65502,16},{65503,16},{65504,16},{65505,16},{65506,16},{65507,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {2041,11},{65508,16},{65509,16},{65510,16},{65511,16},{65512,16},{65513,16},{65514,16},{65515,16},{65516,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, + {16352,14},{65517,16},{65518,16},{65519,16},{65520,16},{65521,16},{65522,16},{65523,16},{65524,16},{65525,16},{0,0},{0,0},{0,0},{0,0},{0,0}, + {1018,10},{32707,15},{65526,16},{65527,16},{65528,16},{65529,16},{65530,16},{65531,16},{65532,16},{65533,16},{65534,16},{0,0},{0,0},{0,0},{0,0},{0,0} + }; + static const int YQT[] = {16,11,10,16,24,40,51,61,12,12,14,19,26,58,60,55,14,13,16,24,40,57,69,56,14,17,22,29,51,87,80,62,18,22, + 37,56,68,109,103,77,24,35,55,64,81,104,113,92,49,64,78,87,103,121,120,101,72,92,95,98,112,100,103,99}; + static const int UVQT[] = {17,18,24,47,99,99,99,99,18,21,26,66,99,99,99,99,24,26,56,99,99,99,99,99,47,66,99,99,99,99,99,99, + 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99}; + static const float aasf[] = { 1.0f * 2.828427125f, 1.387039845f * 2.828427125f, 1.306562965f * 2.828427125f, 1.175875602f * 2.828427125f, + 1.0f * 2.828427125f, 0.785694958f * 2.828427125f, 0.541196100f * 2.828427125f, 0.275899379f * 2.828427125f }; + + int row, col, i, k; + float fdtbl_Y[64], fdtbl_UV[64]; + unsigned char YTable[64], UVTable[64]; + + if(!data || !width || !height || comp > 4 || comp < 1) { + return 0; + } + + quality = quality ? quality : 90; + quality = quality < 1 ? 1 : quality > 100 ? 100 : quality; + quality = quality < 50 ? 5000 / quality : 200 - quality * 2; + + for(i = 0; i < 64; ++i) { + int uvti, yti = (YQT[i]*quality+50)/100; + YTable[stbiw__jpg_ZigZag[i]] = (unsigned char) (yti < 1 ? 1 : yti > 255 ? 255 : yti); + uvti = (UVQT[i]*quality+50)/100; + UVTable[stbiw__jpg_ZigZag[i]] = (unsigned char) (uvti < 1 ? 1 : uvti > 255 ? 255 : uvti); + } + + for(row = 0, k = 0; row < 8; ++row) { + for(col = 0; col < 8; ++col, ++k) { + fdtbl_Y[k] = 1 / (YTable [stbiw__jpg_ZigZag[k]] * aasf[row] * aasf[col]); + fdtbl_UV[k] = 1 / (UVTable[stbiw__jpg_ZigZag[k]] * aasf[row] * aasf[col]); + } + } + + // Write Headers + { + static const unsigned char head0[] = { 0xFF,0xD8,0xFF,0xE0,0,0x10,'J','F','I','F',0,1,1,0,0,1,0,1,0,0,0xFF,0xDB,0,0x84,0 }; + static const unsigned char head2[] = { 0xFF,0xDA,0,0xC,3,1,0,2,0x11,3,0x11,0,0x3F,0 }; + const unsigned char head1[] = { 0xFF,0xC0,0,0x11,8,(unsigned char)(height>>8),STBIW_UCHAR(height),(unsigned char)(width>>8),STBIW_UCHAR(width), + 3,1,0x11,0,2,0x11,1,3,0x11,1,0xFF,0xC4,0x01,0xA2,0 }; + s->func(s->context, (void*)head0, sizeof(head0)); + s->func(s->context, (void*)YTable, sizeof(YTable)); + stbiw__putc(s, 1); + s->func(s->context, UVTable, sizeof(UVTable)); + s->func(s->context, (void*)head1, sizeof(head1)); + s->func(s->context, (void*)(std_dc_luminance_nrcodes+1), sizeof(std_dc_luminance_nrcodes)-1); + s->func(s->context, (void*)std_dc_luminance_values, sizeof(std_dc_luminance_values)); + stbiw__putc(s, 0x10); // HTYACinfo + s->func(s->context, (void*)(std_ac_luminance_nrcodes+1), sizeof(std_ac_luminance_nrcodes)-1); + s->func(s->context, (void*)std_ac_luminance_values, sizeof(std_ac_luminance_values)); + stbiw__putc(s, 1); // HTUDCinfo + s->func(s->context, (void*)(std_dc_chrominance_nrcodes+1), sizeof(std_dc_chrominance_nrcodes)-1); + s->func(s->context, (void*)std_dc_chrominance_values, sizeof(std_dc_chrominance_values)); + stbiw__putc(s, 0x11); // HTUACinfo + s->func(s->context, (void*)(std_ac_chrominance_nrcodes+1), sizeof(std_ac_chrominance_nrcodes)-1); + s->func(s->context, (void*)std_ac_chrominance_values, sizeof(std_ac_chrominance_values)); + s->func(s->context, (void*)head2, sizeof(head2)); + } + + // Encode 8x8 macroblocks + { + static const unsigned short fillBits[] = {0x7F, 7}; + const unsigned char *imageData = (const unsigned char *)data; + int DCY=0, DCU=0, DCV=0; + int bitBuf=0, bitCnt=0; + // comp == 2 is grey+alpha (alpha is ignored) + int ofsG = comp > 2 ? 1 : 0, ofsB = comp > 2 ? 2 : 0; + int x, y, pos; + for(y = 0; y < height; y += 8) { + for(x = 0; x < width; x += 8) { + float YDU[64], UDU[64], VDU[64]; + for(row = y, pos = 0; row < y+8; ++row) { + for(col = x; col < x+8; ++col, ++pos) { + int p = (stbi__flip_vertically_on_write ? height-1-row : row)*width*comp + col*comp; + float r, g, b; + if(row >= height) { + p -= width*comp*(row+1 - height); + } + if(col >= width) { + p -= comp*(col+1 - width); + } + + r = imageData[p+0]; + g = imageData[p+ofsG]; + b = imageData[p+ofsB]; + YDU[pos]=+0.29900f*r+0.58700f*g+0.11400f*b-128; + UDU[pos]=-0.16874f*r-0.33126f*g+0.50000f*b; + VDU[pos]=+0.50000f*r-0.41869f*g-0.08131f*b; + } + } + + DCY = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, YDU, fdtbl_Y, DCY, YDC_HT, YAC_HT); + DCU = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, UDU, fdtbl_UV, DCU, UVDC_HT, UVAC_HT); + DCV = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, VDU, fdtbl_UV, DCV, UVDC_HT, UVAC_HT); + } + } + + // Do the bit alignment of the EOI marker + stbiw__jpg_writeBits(s, &bitBuf, &bitCnt, fillBits); + } + + // EOI + stbiw__putc(s, 0xFF); + stbiw__putc(s, 0xD9); + + return 1; +} + +STBIWDEF int stbi_write_jpg_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data, int quality) +{ + stbi__write_context s; + stbi__start_write_callbacks(&s, func, context); + return stbi_write_jpg_core(&s, x, y, comp, (void *) data, quality); +} + + +#ifndef STBI_WRITE_NO_STDIO +STBIWDEF int stbi_write_jpg(char const *filename, int x, int y, int comp, const void *data, int quality) +{ + stbi__write_context s; + if (stbi__start_write_file(&s,filename)) { + int r = stbi_write_jpg_core(&s, x, y, comp, data, quality); + stbi__end_write_file(&s); + return r; + } else + return 0; +} +#endif + +#endif // STB_IMAGE_WRITE_IMPLEMENTATION + +/* Revision history + 1.09 (2018-02-11) + fix typo in zlib quality API, improve STB_I_W_STATIC in C++ + 1.08 (2018-01-29) + add stbi__flip_vertically_on_write, external zlib, zlib quality, choose PNG filter + 1.07 (2017-07-24) + doc fix + 1.06 (2017-07-23) + writing JPEG (using Jon Olick's code) + 1.05 ??? + 1.04 (2017-03-03) + monochrome BMP expansion + 1.03 ??? + 1.02 (2016-04-02) + avoid allocating large structures on the stack + 1.01 (2016-01-16) + STBIW_REALLOC_SIZED: support allocators with no realloc support + avoid race-condition in crc initialization + minor compile issues + 1.00 (2015-09-14) + installable file IO function + 0.99 (2015-09-13) + warning fixes; TGA rle support + 0.98 (2015-04-08) + added STBIW_MALLOC, STBIW_ASSERT etc + 0.97 (2015-01-18) + fixed HDR asserts, rewrote HDR rle logic + 0.96 (2015-01-17) + add HDR output + fix monochrome BMP + 0.95 (2014-08-17) + add monochrome TGA output + 0.94 (2014-05-31) + rename private functions to avoid conflicts with stb_image.h + 0.93 (2014-05-27) + warning fixes + 0.92 (2010-08-01) + casts to unsigned char to fix warnings + 0.91 (2010-07-17) + first public release + 0.90 first internal release +*/ + +/* +------------------------------------------------------------------------------ +This software is available under 2 licenses -- choose whichever you prefer. +------------------------------------------------------------------------------ +ALTERNATIVE A - MIT License +Copyright (c) 2017 Sean Barrett +Permission is hereby granted, free of charge, to any person obtaining a copy of +this software and associated documentation files (the "Software"), to deal in +the Software without restriction, including without limitation the rights to +use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies +of the Software, and to permit persons to whom the Software is furnished to do +so, subject to the following conditions: +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +------------------------------------------------------------------------------ +ALTERNATIVE B - Public Domain (www.unlicense.org) +This is free and unencumbered software released into the public domain. +Anyone is free to copy, modify, publish, use, compile, sell, or distribute this +software, either in source code form or as a compiled binary, for any purpose, +commercial or non-commercial, and by any means. +In jurisdictions that recognize copyright laws, the author or authors of this +software dedicate any and all copyright interest in the software to the public +domain. We make this dedication for the benefit of the public at large and to +the detriment of our heirs and successors. We intend this dedication to be an +overt act of relinquishment in perpetuity of all present and future rights to +this software under copyright law. +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN +ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION +WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +------------------------------------------------------------------------------ +*/ diff --git a/src/util.h b/src/util.h new file mode 100644 index 0000000..821c3e6 --- /dev/null +++ b/src/util.h @@ -0,0 +1,49 @@ +#ifndef UTIL_H +#define UTIL_H + +#include // fprintf() +#include // exit() + +#define MIN(a, b) ((a) < (b) ? (a) : (b)) +#define MAX(a, b) ((a) > (b) ? (a) : (b)) +#define UNUSED(a) ((void) (a)) + +#define D(...) do { \ + fprintf(stderr, "DEBUG: %s(): ", __func__); \ + fprintf(stderr, __VA_ARGS__); \ + fputs("\n", stderr); \ +} while (0) + +#define die(...) do { \ + fputs("FATAL: ", stderr); \ + fprintf(stderr, __VA_ARGS__); \ + fputs("\n", stderr); \ + exit(EXIT_FAILURE); \ +} while (0) + +// calculate squared euclidean distance between two points +static inline float +distance_squared( + const size_t num_floats, + const float * const a, + const float * const b +) { + float r = 0.0; + + for (size_t i = 0; i < num_floats; i++) { + r += (b[i] - a[i]) * (b[i] - a[i]); + } + + // return squared distance + return r; +} + +static inline float +km_score( + const float mean_distance, + const size_t num_empty +) { + return 1.0 / (mean_distance + num_empty); +} + +#endif /* UTIL_H */ diff --git a/stb_image_write.h b/stb_image_write.h deleted file mode 100644 index c05e958..0000000 --- a/stb_image_write.h +++ /dev/null @@ -1,1568 +0,0 @@ -/* stb_image_write - v1.09 - public domain - http://nothings.org/stb/stb_image_write.h - writes out PNG/BMP/TGA/JPEG/HDR images to C stdio - Sean Barrett 2010-2015 - no warranty implied; use at your own risk - - Before #including, - - #define STB_IMAGE_WRITE_IMPLEMENTATION - - in the file that you want to have the implementation. - - Will probably not work correctly with strict-aliasing optimizations. - - If using a modern Microsoft Compiler, non-safe versions of CRT calls may cause - compilation warnings or even errors. To avoid this, also before #including, - - #define STBI_MSC_SECURE_CRT - -ABOUT: - - This header file is a library for writing images to C stdio. It could be - adapted to write to memory or a general streaming interface; let me know. - - The PNG output is not optimal; it is 20-50% larger than the file - written by a decent optimizing implementation; though providing a custom - zlib compress function (see STBIW_ZLIB_COMPRESS) can mitigate that. - This library is designed for source code compactness and simplicity, - not optimal image file size or run-time performance. - -BUILDING: - - You can #define STBIW_ASSERT(x) before the #include to avoid using assert.h. - You can #define STBIW_MALLOC(), STBIW_REALLOC(), and STBIW_FREE() to replace - malloc,realloc,free. - You can #define STBIW_MEMMOVE() to replace memmove() - You can #define STBIW_ZLIB_COMPRESS to use a custom zlib-style compress function - for PNG compression (instead of the builtin one), it must have the following signature: - unsigned char * my_compress(unsigned char *data, int data_len, int *out_len, int quality); - The returned data will be freed with STBIW_FREE() (free() by default), - so it must be heap allocated with STBIW_MALLOC() (malloc() by default), - -USAGE: - - There are five functions, one for each image file format: - - int stbi_write_png(char const *filename, int w, int h, int comp, const void *data, int stride_in_bytes); - int stbi_write_bmp(char const *filename, int w, int h, int comp, const void *data); - int stbi_write_tga(char const *filename, int w, int h, int comp, const void *data); - int stbi_write_jpg(char const *filename, int w, int h, int comp, const void *data, int quality); - int stbi_write_hdr(char const *filename, int w, int h, int comp, const float *data); - - void stbi_flip_vertically_on_write(int flag); // flag is non-zero to flip data vertically - - There are also five equivalent functions that use an arbitrary write function. You are - expected to open/close your file-equivalent before and after calling these: - - int stbi_write_png_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data, int stride_in_bytes); - int stbi_write_bmp_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data); - int stbi_write_tga_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data); - int stbi_write_hdr_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const float *data); - int stbi_write_jpg_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data, int quality); - - where the callback is: - void stbi_write_func(void *context, void *data, int size); - - You can configure it with these global variables: - int stbi_write_tga_with_rle; // defaults to true; set to 0 to disable RLE - int stbi_write_png_compression_level; // defaults to 8; set to higher for more compression - int stbi_write_force_png_filter; // defaults to -1; set to 0..5 to force a filter mode - - - You can define STBI_WRITE_NO_STDIO to disable the file variant of these - functions, so the library will not use stdio.h at all. However, this will - also disable HDR writing, because it requires stdio for formatted output. - - Each function returns 0 on failure and non-0 on success. - - The functions create an image file defined by the parameters. The image - is a rectangle of pixels stored from left-to-right, top-to-bottom. - Each pixel contains 'comp' channels of data stored interleaved with 8-bits - per channel, in the following order: 1=Y, 2=YA, 3=RGB, 4=RGBA. (Y is - monochrome color.) The rectangle is 'w' pixels wide and 'h' pixels tall. - The *data pointer points to the first byte of the top-left-most pixel. - For PNG, "stride_in_bytes" is the distance in bytes from the first byte of - a row of pixels to the first byte of the next row of pixels. - - PNG creates output files with the same number of components as the input. - The BMP format expands Y to RGB in the file format and does not - output alpha. - - PNG supports writing rectangles of data even when the bytes storing rows of - data are not consecutive in memory (e.g. sub-rectangles of a larger image), - by supplying the stride between the beginning of adjacent rows. The other - formats do not. (Thus you cannot write a native-format BMP through the BMP - writer, both because it is in BGR order and because it may have padding - at the end of the line.) - - PNG allows you to set the deflate compression level by setting the global - variable 'stbi_write_png_compression_level' (it defaults to 8). - - HDR expects linear float data. Since the format is always 32-bit rgb(e) - data, alpha (if provided) is discarded, and for monochrome data it is - replicated across all three channels. - - TGA supports RLE or non-RLE compressed data. To use non-RLE-compressed - data, set the global variable 'stbi_write_tga_with_rle' to 0. - - JPEG does ignore alpha channels in input data; quality is between 1 and 100. - Higher quality looks better but results in a bigger image. - JPEG baseline (no JPEG progressive). - -CREDITS: - - - Sean Barrett - PNG/BMP/TGA - Baldur Karlsson - HDR - Jean-Sebastien Guay - TGA monochrome - Tim Kelsey - misc enhancements - Alan Hickman - TGA RLE - Emmanuel Julien - initial file IO callback implementation - Jon Olick - original jo_jpeg.cpp code - Daniel Gibson - integrate JPEG, allow external zlib - Aarni Koskela - allow choosing PNG filter - - bugfixes: - github:Chribba - Guillaume Chereau - github:jry2 - github:romigrou - Sergio Gonzalez - Jonas Karlsson - Filip Wasil - Thatcher Ulrich - github:poppolopoppo - Patrick Boettcher - github:xeekworx - Cap Petschulat - Simon Rodriguez - Ivan Tikhonov - github:ignotion - Adam Schackart - -LICENSE - - See end of file for license information. - -*/ - -#ifndef INCLUDE_STB_IMAGE_WRITE_H -#define INCLUDE_STB_IMAGE_WRITE_H - -// if STB_IMAGE_WRITE_STATIC causes problems, try defining STBIWDEF to 'inline' or 'static inline' -#ifndef STBIWDEF -#ifdef STB_IMAGE_WRITE_STATIC -#define STBIWDEF static -#else -#ifdef __cplusplus -#define STBIWDEF extern "C" -#else -#define STBIWDEF extern -#endif -#endif -#endif - -#ifndef STB_IMAGE_WRITE_STATIC // C++ forbids static forward declarations -extern int stbi_write_tga_with_rle; -extern int stbi_write_png_compression_level; -extern int stbi_write_force_png_filter; -#endif - -#ifndef STBI_WRITE_NO_STDIO -STBIWDEF int stbi_write_png(char const *filename, int w, int h, int comp, const void *data, int stride_in_bytes); -STBIWDEF int stbi_write_bmp(char const *filename, int w, int h, int comp, const void *data); -STBIWDEF int stbi_write_tga(char const *filename, int w, int h, int comp, const void *data); -STBIWDEF int stbi_write_hdr(char const *filename, int w, int h, int comp, const float *data); -STBIWDEF int stbi_write_jpg(char const *filename, int x, int y, int comp, const void *data, int quality); -#endif - -typedef void stbi_write_func(void *context, void *data, int size); - -STBIWDEF int stbi_write_png_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data, int stride_in_bytes); -STBIWDEF int stbi_write_bmp_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data); -STBIWDEF int stbi_write_tga_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const void *data); -STBIWDEF int stbi_write_hdr_to_func(stbi_write_func *func, void *context, int w, int h, int comp, const float *data); -STBIWDEF int stbi_write_jpg_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data, int quality); - -STBIWDEF void stbi_flip_vertically_on_write(int flip_boolean); - -#endif//INCLUDE_STB_IMAGE_WRITE_H - -#ifdef STB_IMAGE_WRITE_IMPLEMENTATION - -#ifdef _WIN32 - #ifndef _CRT_SECURE_NO_WARNINGS - #define _CRT_SECURE_NO_WARNINGS - #endif - #ifndef _CRT_NONSTDC_NO_DEPRECATE - #define _CRT_NONSTDC_NO_DEPRECATE - #endif -#endif - -#ifndef STBI_WRITE_NO_STDIO -#include -#endif // STBI_WRITE_NO_STDIO - -#include -#include -#include -#include - -#if defined(STBIW_MALLOC) && defined(STBIW_FREE) && (defined(STBIW_REALLOC) || defined(STBIW_REALLOC_SIZED)) -// ok -#elif !defined(STBIW_MALLOC) && !defined(STBIW_FREE) && !defined(STBIW_REALLOC) && !defined(STBIW_REALLOC_SIZED) -// ok -#else -#error "Must define all or none of STBIW_MALLOC, STBIW_FREE, and STBIW_REALLOC (or STBIW_REALLOC_SIZED)." -#endif - -#ifndef STBIW_MALLOC -#define STBIW_MALLOC(sz) malloc(sz) -#define STBIW_REALLOC(p,newsz) realloc(p,newsz) -#define STBIW_FREE(p) free(p) -#endif - -#ifndef STBIW_REALLOC_SIZED -#define STBIW_REALLOC_SIZED(p,oldsz,newsz) STBIW_REALLOC(p,newsz) -#endif - - -#ifndef STBIW_MEMMOVE -#define STBIW_MEMMOVE(a,b,sz) memmove(a,b,sz) -#endif - - -#ifndef STBIW_ASSERT -#include -#define STBIW_ASSERT(x) assert(x) -#endif - -#define STBIW_UCHAR(x) (unsigned char) ((x) & 0xff) - -#ifdef STB_IMAGE_WRITE_STATIC -static int stbi__flip_vertically_on_write=0; -static int stbi_write_png_compression_level = 8; -static int stbi_write_tga_with_rle = 1; -static int stbi_write_force_png_filter = -1; -#else -int stbi_write_png_compression_level = 8; -int stbi__flip_vertically_on_write=0; -int stbi_write_tga_with_rle = 1; -int stbi_write_force_png_filter = -1; -#endif - -STBIWDEF void stbi_flip_vertically_on_write(int flag) -{ - stbi__flip_vertically_on_write = flag; -} - -typedef struct -{ - stbi_write_func *func; - void *context; -} stbi__write_context; - -// initialize a callback-based context -static void stbi__start_write_callbacks(stbi__write_context *s, stbi_write_func *c, void *context) -{ - s->func = c; - s->context = context; -} - -#ifndef STBI_WRITE_NO_STDIO - -static void stbi__stdio_write(void *context, void *data, int size) -{ - fwrite(data,1,size,(FILE*) context); -} - -static int stbi__start_write_file(stbi__write_context *s, const char *filename) -{ - FILE *f; -#ifdef STBI_MSC_SECURE_CRT - if (fopen_s(&f, filename, "wb")) - f = NULL; -#else - f = fopen(filename, "wb"); -#endif - stbi__start_write_callbacks(s, stbi__stdio_write, (void *) f); - return f != NULL; -} - -static void stbi__end_write_file(stbi__write_context *s) -{ - fclose((FILE *)s->context); -} - -#endif // !STBI_WRITE_NO_STDIO - -typedef unsigned int stbiw_uint32; -typedef int stb_image_write_test[sizeof(stbiw_uint32)==4 ? 1 : -1]; - -static void stbiw__writefv(stbi__write_context *s, const char *fmt, va_list v) -{ - while (*fmt) { - switch (*fmt++) { - case ' ': break; - case '1': { unsigned char x = STBIW_UCHAR(va_arg(v, int)); - s->func(s->context,&x,1); - break; } - case '2': { int x = va_arg(v,int); - unsigned char b[2]; - b[0] = STBIW_UCHAR(x); - b[1] = STBIW_UCHAR(x>>8); - s->func(s->context,b,2); - break; } - case '4': { stbiw_uint32 x = va_arg(v,int); - unsigned char b[4]; - b[0]=STBIW_UCHAR(x); - b[1]=STBIW_UCHAR(x>>8); - b[2]=STBIW_UCHAR(x>>16); - b[3]=STBIW_UCHAR(x>>24); - s->func(s->context,b,4); - break; } - default: - STBIW_ASSERT(0); - return; - } - } -} - -static void stbiw__writef(stbi__write_context *s, const char *fmt, ...) -{ - va_list v; - va_start(v, fmt); - stbiw__writefv(s, fmt, v); - va_end(v); -} - -static void stbiw__putc(stbi__write_context *s, unsigned char c) -{ - s->func(s->context, &c, 1); -} - -static void stbiw__write3(stbi__write_context *s, unsigned char a, unsigned char b, unsigned char c) -{ - unsigned char arr[3]; - arr[0] = a, arr[1] = b, arr[2] = c; - s->func(s->context, arr, 3); -} - -static void stbiw__write_pixel(stbi__write_context *s, int rgb_dir, int comp, int write_alpha, int expand_mono, unsigned char *d) -{ - unsigned char bg[3] = { 255, 0, 255}, px[3]; - int k; - - if (write_alpha < 0) - s->func(s->context, &d[comp - 1], 1); - - switch (comp) { - case 2: // 2 pixels = mono + alpha, alpha is written separately, so same as 1-channel case - case 1: - if (expand_mono) - stbiw__write3(s, d[0], d[0], d[0]); // monochrome bmp - else - s->func(s->context, d, 1); // monochrome TGA - break; - case 4: - if (!write_alpha) { - // composite against pink background - for (k = 0; k < 3; ++k) - px[k] = bg[k] + ((d[k] - bg[k]) * d[3]) / 255; - stbiw__write3(s, px[1 - rgb_dir], px[1], px[1 + rgb_dir]); - break; - } - /* FALLTHROUGH */ - case 3: - stbiw__write3(s, d[1 - rgb_dir], d[1], d[1 + rgb_dir]); - break; - } - if (write_alpha > 0) - s->func(s->context, &d[comp - 1], 1); -} - -static void stbiw__write_pixels(stbi__write_context *s, int rgb_dir, int vdir, int x, int y, int comp, void *data, int write_alpha, int scanline_pad, int expand_mono) -{ - stbiw_uint32 zero = 0; - int i,j, j_end; - - if (y <= 0) - return; - - if (stbi__flip_vertically_on_write) - vdir *= -1; - - if (vdir < 0) - j_end = -1, j = y-1; - else - j_end = y, j = 0; - - for (; j != j_end; j += vdir) { - for (i=0; i < x; ++i) { - unsigned char *d = (unsigned char *) data + (j*x+i)*comp; - stbiw__write_pixel(s, rgb_dir, comp, write_alpha, expand_mono, d); - } - s->func(s->context, &zero, scanline_pad); - } -} - -static int stbiw__outfile(stbi__write_context *s, int rgb_dir, int vdir, int x, int y, int comp, int expand_mono, void *data, int alpha, int pad, const char *fmt, ...) -{ - if (y < 0 || x < 0) { - return 0; - } else { - va_list v; - va_start(v, fmt); - stbiw__writefv(s, fmt, v); - va_end(v); - stbiw__write_pixels(s,rgb_dir,vdir,x,y,comp,data,alpha,pad, expand_mono); - return 1; - } -} - -static int stbi_write_bmp_core(stbi__write_context *s, int x, int y, int comp, const void *data) -{ - int pad = (-x*3) & 3; - return stbiw__outfile(s,-1,-1,x,y,comp,1,(void *) data,0,pad, - "11 4 22 4" "4 44 22 444444", - 'B', 'M', 14+40+(x*3+pad)*y, 0,0, 14+40, // file header - 40, x,y, 1,24, 0,0,0,0,0,0); // bitmap header -} - -STBIWDEF int stbi_write_bmp_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data) -{ - stbi__write_context s; - stbi__start_write_callbacks(&s, func, context); - return stbi_write_bmp_core(&s, x, y, comp, data); -} - -#ifndef STBI_WRITE_NO_STDIO -STBIWDEF int stbi_write_bmp(char const *filename, int x, int y, int comp, const void *data) -{ - stbi__write_context s; - if (stbi__start_write_file(&s,filename)) { - int r = stbi_write_bmp_core(&s, x, y, comp, data); - stbi__end_write_file(&s); - return r; - } else - return 0; -} -#endif //!STBI_WRITE_NO_STDIO - -static int stbi_write_tga_core(stbi__write_context *s, int x, int y, int comp, void *data) -{ - int has_alpha = (comp == 2 || comp == 4); - int colorbytes = has_alpha ? comp-1 : comp; - int format = colorbytes < 2 ? 3 : 2; // 3 color channels (RGB/RGBA) = 2, 1 color channel (Y/YA) = 3 - - if (y < 0 || x < 0) - return 0; - - if (!stbi_write_tga_with_rle) { - return stbiw__outfile(s, -1, -1, x, y, comp, 0, (void *) data, has_alpha, 0, - "111 221 2222 11", 0, 0, format, 0, 0, 0, 0, 0, x, y, (colorbytes + has_alpha) * 8, has_alpha * 8); - } else { - int i,j,k; - int jend, jdir; - - stbiw__writef(s, "111 221 2222 11", 0,0,format+8, 0,0,0, 0,0,x,y, (colorbytes + has_alpha) * 8, has_alpha * 8); - - if (stbi__flip_vertically_on_write) { - j = 0; - jend = y; - jdir = 1; - } else { - j = y-1; - jend = -1; - jdir = -1; - } - for (; j != jend; j += jdir) { - unsigned char *row = (unsigned char *) data + j * x * comp; - int len; - - for (i = 0; i < x; i += len) { - unsigned char *begin = row + i * comp; - int diff = 1; - len = 1; - - if (i < x - 1) { - ++len; - diff = memcmp(begin, row + (i + 1) * comp, comp); - if (diff) { - const unsigned char *prev = begin; - for (k = i + 2; k < x && len < 128; ++k) { - if (memcmp(prev, row + k * comp, comp)) { - prev += comp; - ++len; - } else { - --len; - break; - } - } - } else { - for (k = i + 2; k < x && len < 128; ++k) { - if (!memcmp(begin, row + k * comp, comp)) { - ++len; - } else { - break; - } - } - } - } - - if (diff) { - unsigned char header = STBIW_UCHAR(len - 1); - s->func(s->context, &header, 1); - for (k = 0; k < len; ++k) { - stbiw__write_pixel(s, -1, comp, has_alpha, 0, begin + k * comp); - } - } else { - unsigned char header = STBIW_UCHAR(len - 129); - s->func(s->context, &header, 1); - stbiw__write_pixel(s, -1, comp, has_alpha, 0, begin); - } - } - } - } - return 1; -} - -STBIWDEF int stbi_write_tga_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data) -{ - stbi__write_context s; - stbi__start_write_callbacks(&s, func, context); - return stbi_write_tga_core(&s, x, y, comp, (void *) data); -} - -#ifndef STBI_WRITE_NO_STDIO -STBIWDEF int stbi_write_tga(char const *filename, int x, int y, int comp, const void *data) -{ - stbi__write_context s; - if (stbi__start_write_file(&s,filename)) { - int r = stbi_write_tga_core(&s, x, y, comp, (void *) data); - stbi__end_write_file(&s); - return r; - } else - return 0; -} -#endif - -// ************************************************************************************************* -// Radiance RGBE HDR writer -// by Baldur Karlsson - -#define stbiw__max(a, b) ((a) > (b) ? (a) : (b)) - -void stbiw__linear_to_rgbe(unsigned char *rgbe, float *linear) -{ - int exponent; - float maxcomp = stbiw__max(linear[0], stbiw__max(linear[1], linear[2])); - - if (maxcomp < 1e-32f) { - rgbe[0] = rgbe[1] = rgbe[2] = rgbe[3] = 0; - } else { - float normalize = (float) frexp(maxcomp, &exponent) * 256.0f/maxcomp; - - rgbe[0] = (unsigned char)(linear[0] * normalize); - rgbe[1] = (unsigned char)(linear[1] * normalize); - rgbe[2] = (unsigned char)(linear[2] * normalize); - rgbe[3] = (unsigned char)(exponent + 128); - } -} - -void stbiw__write_run_data(stbi__write_context *s, int length, unsigned char databyte) -{ - unsigned char lengthbyte = STBIW_UCHAR(length+128); - STBIW_ASSERT(length+128 <= 255); - s->func(s->context, &lengthbyte, 1); - s->func(s->context, &databyte, 1); -} - -void stbiw__write_dump_data(stbi__write_context *s, int length, unsigned char *data) -{ - unsigned char lengthbyte = STBIW_UCHAR(length); - STBIW_ASSERT(length <= 128); // inconsistent with spec but consistent with official code - s->func(s->context, &lengthbyte, 1); - s->func(s->context, data, length); -} - -void stbiw__write_hdr_scanline(stbi__write_context *s, int width, int ncomp, unsigned char *scratch, float *scanline) -{ - unsigned char scanlineheader[4] = { 2, 2, 0, 0 }; - unsigned char rgbe[4]; - float linear[3]; - int x; - - scanlineheader[2] = (width&0xff00)>>8; - scanlineheader[3] = (width&0x00ff); - - /* skip RLE for images too small or large */ - if (width < 8 || width >= 32768) { - for (x=0; x < width; x++) { - switch (ncomp) { - case 4: /* fallthrough */ - case 3: linear[2] = scanline[x*ncomp + 2]; - linear[1] = scanline[x*ncomp + 1]; - linear[0] = scanline[x*ncomp + 0]; - break; - default: - linear[0] = linear[1] = linear[2] = scanline[x*ncomp + 0]; - break; - } - stbiw__linear_to_rgbe(rgbe, linear); - s->func(s->context, rgbe, 4); - } - } else { - int c,r; - /* encode into scratch buffer */ - for (x=0; x < width; x++) { - switch(ncomp) { - case 4: /* fallthrough */ - case 3: linear[2] = scanline[x*ncomp + 2]; - linear[1] = scanline[x*ncomp + 1]; - linear[0] = scanline[x*ncomp + 0]; - break; - default: - linear[0] = linear[1] = linear[2] = scanline[x*ncomp + 0]; - break; - } - stbiw__linear_to_rgbe(rgbe, linear); - scratch[x + width*0] = rgbe[0]; - scratch[x + width*1] = rgbe[1]; - scratch[x + width*2] = rgbe[2]; - scratch[x + width*3] = rgbe[3]; - } - - s->func(s->context, scanlineheader, 4); - - /* RLE each component separately */ - for (c=0; c < 4; c++) { - unsigned char *comp = &scratch[width*c]; - - x = 0; - while (x < width) { - // find first run - r = x; - while (r+2 < width) { - if (comp[r] == comp[r+1] && comp[r] == comp[r+2]) - break; - ++r; - } - if (r+2 >= width) - r = width; - // dump up to first run - while (x < r) { - int len = r-x; - if (len > 128) len = 128; - stbiw__write_dump_data(s, len, &comp[x]); - x += len; - } - // if there's a run, output it - if (r+2 < width) { // same test as what we break out of in search loop, so only true if we break'd - // find next byte after run - while (r < width && comp[r] == comp[x]) - ++r; - // output run up to r - while (x < r) { - int len = r-x; - if (len > 127) len = 127; - stbiw__write_run_data(s, len, comp[x]); - x += len; - } - } - } - } - } -} - -static int stbi_write_hdr_core(stbi__write_context *s, int x, int y, int comp, float *data) -{ - if (y <= 0 || x <= 0 || data == NULL) - return 0; - else { - // Each component is stored separately. Allocate scratch space for full output scanline. - unsigned char *scratch = (unsigned char *) STBIW_MALLOC(x*4); - int i, len; - char buffer[128]; - char header[] = "#?RADIANCE\n# Written by stb_image_write.h\nFORMAT=32-bit_rle_rgbe\n"; - s->func(s->context, header, sizeof(header)-1); - -#ifdef STBI_MSC_SECURE_CRT - len = sprintf_s(buffer, "EXPOSURE= 1.0000000000000\n\n-Y %d +X %d\n", y, x); -#else - len = sprintf(buffer, "EXPOSURE= 1.0000000000000\n\n-Y %d +X %d\n", y, x); -#endif - s->func(s->context, buffer, len); - - for(i=0; i < y; i++) - stbiw__write_hdr_scanline(s, x, comp, scratch, data + comp*x*(stbi__flip_vertically_on_write ? y-1-i : i)*x); - STBIW_FREE(scratch); - return 1; - } -} - -STBIWDEF int stbi_write_hdr_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const float *data) -{ - stbi__write_context s; - stbi__start_write_callbacks(&s, func, context); - return stbi_write_hdr_core(&s, x, y, comp, (float *) data); -} - -#ifndef STBI_WRITE_NO_STDIO -STBIWDEF int stbi_write_hdr(char const *filename, int x, int y, int comp, const float *data) -{ - stbi__write_context s; - if (stbi__start_write_file(&s,filename)) { - int r = stbi_write_hdr_core(&s, x, y, comp, (float *) data); - stbi__end_write_file(&s); - return r; - } else - return 0; -} -#endif // STBI_WRITE_NO_STDIO - - -////////////////////////////////////////////////////////////////////////////// -// -// PNG writer -// - -#ifndef STBIW_ZLIB_COMPRESS -// stretchy buffer; stbiw__sbpush() == vector<>::push_back() -- stbiw__sbcount() == vector<>::size() -#define stbiw__sbraw(a) ((int *) (a) - 2) -#define stbiw__sbm(a) stbiw__sbraw(a)[0] -#define stbiw__sbn(a) stbiw__sbraw(a)[1] - -#define stbiw__sbneedgrow(a,n) ((a)==0 || stbiw__sbn(a)+n >= stbiw__sbm(a)) -#define stbiw__sbmaybegrow(a,n) (stbiw__sbneedgrow(a,(n)) ? stbiw__sbgrow(a,n) : 0) -#define stbiw__sbgrow(a,n) stbiw__sbgrowf((void **) &(a), (n), sizeof(*(a))) - -#define stbiw__sbpush(a, v) (stbiw__sbmaybegrow(a,1), (a)[stbiw__sbn(a)++] = (v)) -#define stbiw__sbcount(a) ((a) ? stbiw__sbn(a) : 0) -#define stbiw__sbfree(a) ((a) ? STBIW_FREE(stbiw__sbraw(a)),0 : 0) - -static void *stbiw__sbgrowf(void **arr, int increment, int itemsize) -{ - int m = *arr ? 2*stbiw__sbm(*arr)+increment : increment+1; - void *p = STBIW_REALLOC_SIZED(*arr ? stbiw__sbraw(*arr) : 0, *arr ? (stbiw__sbm(*arr)*itemsize + sizeof(int)*2) : 0, itemsize * m + sizeof(int)*2); - STBIW_ASSERT(p); - if (p) { - if (!*arr) ((int *) p)[1] = 0; - *arr = (void *) ((int *) p + 2); - stbiw__sbm(*arr) = m; - } - return *arr; -} - -static unsigned char *stbiw__zlib_flushf(unsigned char *data, unsigned int *bitbuffer, int *bitcount) -{ - while (*bitcount >= 8) { - stbiw__sbpush(data, STBIW_UCHAR(*bitbuffer)); - *bitbuffer >>= 8; - *bitcount -= 8; - } - return data; -} - -static int stbiw__zlib_bitrev(int code, int codebits) -{ - int res=0; - while (codebits--) { - res = (res << 1) | (code & 1); - code >>= 1; - } - return res; -} - -static unsigned int stbiw__zlib_countm(unsigned char *a, unsigned char *b, int limit) -{ - int i; - for (i=0; i < limit && i < 258; ++i) - if (a[i] != b[i]) break; - return i; -} - -static unsigned int stbiw__zhash(unsigned char *data) -{ - stbiw_uint32 hash = data[0] + (data[1] << 8) + (data[2] << 16); - hash ^= hash << 3; - hash += hash >> 5; - hash ^= hash << 4; - hash += hash >> 17; - hash ^= hash << 25; - hash += hash >> 6; - return hash; -} - -#define stbiw__zlib_flush() (out = stbiw__zlib_flushf(out, &bitbuf, &bitcount)) -#define stbiw__zlib_add(code,codebits) \ - (bitbuf |= (code) << bitcount, bitcount += (codebits), stbiw__zlib_flush()) -#define stbiw__zlib_huffa(b,c) stbiw__zlib_add(stbiw__zlib_bitrev(b,c),c) -// default huffman tables -#define stbiw__zlib_huff1(n) stbiw__zlib_huffa(0x30 + (n), 8) -#define stbiw__zlib_huff2(n) stbiw__zlib_huffa(0x190 + (n)-144, 9) -#define stbiw__zlib_huff3(n) stbiw__zlib_huffa(0 + (n)-256,7) -#define stbiw__zlib_huff4(n) stbiw__zlib_huffa(0xc0 + (n)-280,8) -#define stbiw__zlib_huff(n) ((n) <= 143 ? stbiw__zlib_huff1(n) : (n) <= 255 ? stbiw__zlib_huff2(n) : (n) <= 279 ? stbiw__zlib_huff3(n) : stbiw__zlib_huff4(n)) -#define stbiw__zlib_huffb(n) ((n) <= 143 ? stbiw__zlib_huff1(n) : stbiw__zlib_huff2(n)) - -#define stbiw__ZHASH 16384 - -#endif // STBIW_ZLIB_COMPRESS - -unsigned char * stbi_zlib_compress(unsigned char *data, int data_len, int *out_len, int quality) -{ -#ifdef STBIW_ZLIB_COMPRESS - // user provided a zlib compress implementation, use that - return STBIW_ZLIB_COMPRESS(data, data_len, out_len, quality); -#else // use builtin - static unsigned short lengthc[] = { 3,4,5,6,7,8,9,10,11,13,15,17,19,23,27,31,35,43,51,59,67,83,99,115,131,163,195,227,258, 259 }; - static unsigned char lengtheb[]= { 0,0,0,0,0,0,0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 0 }; - static unsigned short distc[] = { 1,2,3,4,5,7,9,13,17,25,33,49,65,97,129,193,257,385,513,769,1025,1537,2049,3073,4097,6145,8193,12289,16385,24577, 32768 }; - static unsigned char disteb[] = { 0,0,0,0,1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,11,11,12,12,13,13 }; - unsigned int bitbuf=0; - int i,j, bitcount=0; - unsigned char *out = NULL; - unsigned char ***hash_table = (unsigned char***) STBIW_MALLOC(stbiw__ZHASH * sizeof(char**)); - if (hash_table == NULL) - return NULL; - if (quality < 5) quality = 5; - - stbiw__sbpush(out, 0x78); // DEFLATE 32K window - stbiw__sbpush(out, 0x5e); // FLEVEL = 1 - stbiw__zlib_add(1,1); // BFINAL = 1 - stbiw__zlib_add(1,2); // BTYPE = 1 -- fixed huffman - - for (i=0; i < stbiw__ZHASH; ++i) - hash_table[i] = NULL; - - i=0; - while (i < data_len-3) { - // hash next 3 bytes of data to be compressed - int h = stbiw__zhash(data+i)&(stbiw__ZHASH-1), best=3; - unsigned char *bestloc = 0; - unsigned char **hlist = hash_table[h]; - int n = stbiw__sbcount(hlist); - for (j=0; j < n; ++j) { - if (hlist[j]-data > i-32768) { // if entry lies within window - int d = stbiw__zlib_countm(hlist[j], data+i, data_len-i); - if (d >= best) best=d,bestloc=hlist[j]; - } - } - // when hash table entry is too long, delete half the entries - if (hash_table[h] && stbiw__sbn(hash_table[h]) == 2*quality) { - STBIW_MEMMOVE(hash_table[h], hash_table[h]+quality, sizeof(hash_table[h][0])*quality); - stbiw__sbn(hash_table[h]) = quality; - } - stbiw__sbpush(hash_table[h],data+i); - - if (bestloc) { - // "lazy matching" - check match at *next* byte, and if it's better, do cur byte as literal - h = stbiw__zhash(data+i+1)&(stbiw__ZHASH-1); - hlist = hash_table[h]; - n = stbiw__sbcount(hlist); - for (j=0; j < n; ++j) { - if (hlist[j]-data > i-32767) { - int e = stbiw__zlib_countm(hlist[j], data+i+1, data_len-i-1); - if (e > best) { // if next match is better, bail on current match - bestloc = NULL; - break; - } - } - } - } - - if (bestloc) { - int d = (int) (data+i - bestloc); // distance back - STBIW_ASSERT(d <= 32767 && best <= 258); - for (j=0; best > lengthc[j+1]-1; ++j); - stbiw__zlib_huff(j+257); - if (lengtheb[j]) stbiw__zlib_add(best - lengthc[j], lengtheb[j]); - for (j=0; d > distc[j+1]-1; ++j); - stbiw__zlib_add(stbiw__zlib_bitrev(j,5),5); - if (disteb[j]) stbiw__zlib_add(d - distc[j], disteb[j]); - i += best; - } else { - stbiw__zlib_huffb(data[i]); - ++i; - } - } - // write out final bytes - for (;i < data_len; ++i) - stbiw__zlib_huffb(data[i]); - stbiw__zlib_huff(256); // end of block - // pad with 0 bits to byte boundary - while (bitcount) - stbiw__zlib_add(0,1); - - for (i=0; i < stbiw__ZHASH; ++i) - (void) stbiw__sbfree(hash_table[i]); - STBIW_FREE(hash_table); - - { - // compute adler32 on input - unsigned int s1=1, s2=0; - int blocklen = (int) (data_len % 5552); - j=0; - while (j < data_len) { - for (i=0; i < blocklen; ++i) s1 += data[j+i], s2 += s1; - s1 %= 65521, s2 %= 65521; - j += blocklen; - blocklen = 5552; - } - stbiw__sbpush(out, STBIW_UCHAR(s2 >> 8)); - stbiw__sbpush(out, STBIW_UCHAR(s2)); - stbiw__sbpush(out, STBIW_UCHAR(s1 >> 8)); - stbiw__sbpush(out, STBIW_UCHAR(s1)); - } - *out_len = stbiw__sbn(out); - // make returned pointer freeable - STBIW_MEMMOVE(stbiw__sbraw(out), out, *out_len); - return (unsigned char *) stbiw__sbraw(out); -#endif // STBIW_ZLIB_COMPRESS -} - -static unsigned int stbiw__crc32(unsigned char *buffer, int len) -{ - static unsigned int crc_table[256] = - { - 0x00000000, 0x77073096, 0xEE0E612C, 0x990951BA, 0x076DC419, 0x706AF48F, 0xE963A535, 0x9E6495A3, - 0x0eDB8832, 0x79DCB8A4, 0xE0D5E91E, 0x97D2D988, 0x09B64C2B, 0x7EB17CBD, 0xE7B82D07, 0x90BF1D91, - 0x1DB71064, 0x6AB020F2, 0xF3B97148, 0x84BE41DE, 0x1ADAD47D, 0x6DDDE4EB, 0xF4D4B551, 0x83D385C7, - 0x136C9856, 0x646BA8C0, 0xFD62F97A, 0x8A65C9EC, 0x14015C4F, 0x63066CD9, 0xFA0F3D63, 0x8D080DF5, - 0x3B6E20C8, 0x4C69105E, 0xD56041E4, 0xA2677172, 0x3C03E4D1, 0x4B04D447, 0xD20D85FD, 0xA50AB56B, - 0x35B5A8FA, 0x42B2986C, 0xDBBBC9D6, 0xACBCF940, 0x32D86CE3, 0x45DF5C75, 0xDCD60DCF, 0xABD13D59, - 0x26D930AC, 0x51DE003A, 0xC8D75180, 0xBFD06116, 0x21B4F4B5, 0x56B3C423, 0xCFBA9599, 0xB8BDA50F, - 0x2802B89E, 0x5F058808, 0xC60CD9B2, 0xB10BE924, 0x2F6F7C87, 0x58684C11, 0xC1611DAB, 0xB6662D3D, - 0x76DC4190, 0x01DB7106, 0x98D220BC, 0xEFD5102A, 0x71B18589, 0x06B6B51F, 0x9FBFE4A5, 0xE8B8D433, - 0x7807C9A2, 0x0F00F934, 0x9609A88E, 0xE10E9818, 0x7F6A0DBB, 0x086D3D2D, 0x91646C97, 0xE6635C01, - 0x6B6B51F4, 0x1C6C6162, 0x856530D8, 0xF262004E, 0x6C0695ED, 0x1B01A57B, 0x8208F4C1, 0xF50FC457, - 0x65B0D9C6, 0x12B7E950, 0x8BBEB8EA, 0xFCB9887C, 0x62DD1DDF, 0x15DA2D49, 0x8CD37CF3, 0xFBD44C65, - 0x4DB26158, 0x3AB551CE, 0xA3BC0074, 0xD4BB30E2, 0x4ADFA541, 0x3DD895D7, 0xA4D1C46D, 0xD3D6F4FB, - 0x4369E96A, 0x346ED9FC, 0xAD678846, 0xDA60B8D0, 0x44042D73, 0x33031DE5, 0xAA0A4C5F, 0xDD0D7CC9, - 0x5005713C, 0x270241AA, 0xBE0B1010, 0xC90C2086, 0x5768B525, 0x206F85B3, 0xB966D409, 0xCE61E49F, - 0x5EDEF90E, 0x29D9C998, 0xB0D09822, 0xC7D7A8B4, 0x59B33D17, 0x2EB40D81, 0xB7BD5C3B, 0xC0BA6CAD, - 0xEDB88320, 0x9ABFB3B6, 0x03B6E20C, 0x74B1D29A, 0xEAD54739, 0x9DD277AF, 0x04DB2615, 0x73DC1683, - 0xE3630B12, 0x94643B84, 0x0D6D6A3E, 0x7A6A5AA8, 0xE40ECF0B, 0x9309FF9D, 0x0A00AE27, 0x7D079EB1, - 0xF00F9344, 0x8708A3D2, 0x1E01F268, 0x6906C2FE, 0xF762575D, 0x806567CB, 0x196C3671, 0x6E6B06E7, - 0xFED41B76, 0x89D32BE0, 0x10DA7A5A, 0x67DD4ACC, 0xF9B9DF6F, 0x8EBEEFF9, 0x17B7BE43, 0x60B08ED5, - 0xD6D6A3E8, 0xA1D1937E, 0x38D8C2C4, 0x4FDFF252, 0xD1BB67F1, 0xA6BC5767, 0x3FB506DD, 0x48B2364B, - 0xD80D2BDA, 0xAF0A1B4C, 0x36034AF6, 0x41047A60, 0xDF60EFC3, 0xA867DF55, 0x316E8EEF, 0x4669BE79, - 0xCB61B38C, 0xBC66831A, 0x256FD2A0, 0x5268E236, 0xCC0C7795, 0xBB0B4703, 0x220216B9, 0x5505262F, - 0xC5BA3BBE, 0xB2BD0B28, 0x2BB45A92, 0x5CB36A04, 0xC2D7FFA7, 0xB5D0CF31, 0x2CD99E8B, 0x5BDEAE1D, - 0x9B64C2B0, 0xEC63F226, 0x756AA39C, 0x026D930A, 0x9C0906A9, 0xEB0E363F, 0x72076785, 0x05005713, - 0x95BF4A82, 0xE2B87A14, 0x7BB12BAE, 0x0CB61B38, 0x92D28E9B, 0xE5D5BE0D, 0x7CDCEFB7, 0x0BDBDF21, - 0x86D3D2D4, 0xF1D4E242, 0x68DDB3F8, 0x1FDA836E, 0x81BE16CD, 0xF6B9265B, 0x6FB077E1, 0x18B74777, - 0x88085AE6, 0xFF0F6A70, 0x66063BCA, 0x11010B5C, 0x8F659EFF, 0xF862AE69, 0x616BFFD3, 0x166CCF45, - 0xA00AE278, 0xD70DD2EE, 0x4E048354, 0x3903B3C2, 0xA7672661, 0xD06016F7, 0x4969474D, 0x3E6E77DB, - 0xAED16A4A, 0xD9D65ADC, 0x40DF0B66, 0x37D83BF0, 0xA9BCAE53, 0xDEBB9EC5, 0x47B2CF7F, 0x30B5FFE9, - 0xBDBDF21C, 0xCABAC28A, 0x53B39330, 0x24B4A3A6, 0xBAD03605, 0xCDD70693, 0x54DE5729, 0x23D967BF, - 0xB3667A2E, 0xC4614AB8, 0x5D681B02, 0x2A6F2B94, 0xB40BBE37, 0xC30C8EA1, 0x5A05DF1B, 0x2D02EF8D - }; - - unsigned int crc = ~0u; - int i; - for (i=0; i < len; ++i) - crc = (crc >> 8) ^ crc_table[buffer[i] ^ (crc & 0xff)]; - return ~crc; -} - -#define stbiw__wpng4(o,a,b,c,d) ((o)[0]=STBIW_UCHAR(a),(o)[1]=STBIW_UCHAR(b),(o)[2]=STBIW_UCHAR(c),(o)[3]=STBIW_UCHAR(d),(o)+=4) -#define stbiw__wp32(data,v) stbiw__wpng4(data, (v)>>24,(v)>>16,(v)>>8,(v)); -#define stbiw__wptag(data,s) stbiw__wpng4(data, s[0],s[1],s[2],s[3]) - -static void stbiw__wpcrc(unsigned char **data, int len) -{ - unsigned int crc = stbiw__crc32(*data - len - 4, len+4); - stbiw__wp32(*data, crc); -} - -static unsigned char stbiw__paeth(int a, int b, int c) -{ - int p = a + b - c, pa = abs(p-a), pb = abs(p-b), pc = abs(p-c); - if (pa <= pb && pa <= pc) return STBIW_UCHAR(a); - if (pb <= pc) return STBIW_UCHAR(b); - return STBIW_UCHAR(c); -} - -// @OPTIMIZE: provide an option that always forces left-predict or paeth predict -static void stbiw__encode_png_line(unsigned char *pixels, int stride_bytes, int width, int height, int y, int n, int filter_type, signed char *line_buffer) -{ - static int mapping[] = { 0,1,2,3,4 }; - static int firstmap[] = { 0,1,0,5,6 }; - int *mymap = (y != 0) ? mapping : firstmap; - int i; - int type = mymap[filter_type]; - unsigned char *z = pixels + stride_bytes * (stbi__flip_vertically_on_write ? height-1-y : y); - int signed_stride = stbi__flip_vertically_on_write ? -stride_bytes : stride_bytes; - for (i = 0; i < n; ++i) { - switch (type) { - case 0: line_buffer[i] = z[i]; break; - case 1: line_buffer[i] = z[i]; break; - case 2: line_buffer[i] = z[i] - z[i-signed_stride]; break; - case 3: line_buffer[i] = z[i] - (z[i-signed_stride]>>1); break; - case 4: line_buffer[i] = (signed char) (z[i] - stbiw__paeth(0,z[i-signed_stride],0)); break; - case 5: line_buffer[i] = z[i]; break; - case 6: line_buffer[i] = z[i]; break; - } - } - for (i=n; i < width*n; ++i) { - switch (type) { - case 0: line_buffer[i] = z[i]; break; - case 1: line_buffer[i] = z[i] - z[i-n]; break; - case 2: line_buffer[i] = z[i] - z[i-signed_stride]; break; - case 3: line_buffer[i] = z[i] - ((z[i-n] + z[i-signed_stride])>>1); break; - case 4: line_buffer[i] = z[i] - stbiw__paeth(z[i-n], z[i-signed_stride], z[i-signed_stride-n]); break; - case 5: line_buffer[i] = z[i] - (z[i-n]>>1); break; - case 6: line_buffer[i] = z[i] - stbiw__paeth(z[i-n], 0,0); break; - } - } -} - -unsigned char *stbi_write_png_to_mem(unsigned char *pixels, int stride_bytes, int x, int y, int n, int *out_len) -{ - int force_filter = stbi_write_force_png_filter; - int ctype[5] = { -1, 0, 4, 2, 6 }; - unsigned char sig[8] = { 137,80,78,71,13,10,26,10 }; - unsigned char *out,*o, *filt, *zlib; - signed char *line_buffer; - int j,zlen; - - if (stride_bytes == 0) - stride_bytes = x * n; - - if (force_filter >= 5) { - force_filter = -1; - } - - filt = (unsigned char *) STBIW_MALLOC((x*n+1) * y); if (!filt) return 0; - line_buffer = (signed char *) STBIW_MALLOC(x * n); if (!line_buffer) { STBIW_FREE(filt); return 0; } - for (j=0; j < y; ++j) { - int filter_type; - if (force_filter > -1) { - filter_type = force_filter; - stbiw__encode_png_line(pixels, stride_bytes, x, y, j, n, force_filter, line_buffer); - } else { // Estimate the best filter by running through all of them: - int best_filter = 0, best_filter_val = 0x7fffffff, est, i; - for (filter_type = 0; filter_type < 5; filter_type++) { - stbiw__encode_png_line(pixels, stride_bytes, x, y, j, n, filter_type, line_buffer); - - // Estimate the entropy of the line using this filter; the less, the better. - est = 0; - for (i = 0; i < x*n; ++i) { - est += abs((signed char) line_buffer[i]); - } - if (est < best_filter_val) { - best_filter_val = est; - best_filter = filter_type; - } - } - if (filter_type != best_filter) { // If the last iteration already got us the best filter, don't redo it - stbiw__encode_png_line(pixels, stride_bytes, x, y, j, n, best_filter, line_buffer); - filter_type = best_filter; - } - } - // when we get here, filter_type contains the filter type, and line_buffer contains the data - filt[j*(x*n+1)] = (unsigned char) filter_type; - STBIW_MEMMOVE(filt+j*(x*n+1)+1, line_buffer, x*n); - } - STBIW_FREE(line_buffer); - zlib = stbi_zlib_compress(filt, y*( x*n+1), &zlen, stbi_write_png_compression_level); - STBIW_FREE(filt); - if (!zlib) return 0; - - // each tag requires 12 bytes of overhead - out = (unsigned char *) STBIW_MALLOC(8 + 12+13 + 12+zlen + 12); - if (!out) return 0; - *out_len = 8 + 12+13 + 12+zlen + 12; - - o=out; - STBIW_MEMMOVE(o,sig,8); o+= 8; - stbiw__wp32(o, 13); // header length - stbiw__wptag(o, "IHDR"); - stbiw__wp32(o, x); - stbiw__wp32(o, y); - *o++ = 8; - *o++ = STBIW_UCHAR(ctype[n]); - *o++ = 0; - *o++ = 0; - *o++ = 0; - stbiw__wpcrc(&o,13); - - stbiw__wp32(o, zlen); - stbiw__wptag(o, "IDAT"); - STBIW_MEMMOVE(o, zlib, zlen); - o += zlen; - STBIW_FREE(zlib); - stbiw__wpcrc(&o, zlen); - - stbiw__wp32(o,0); - stbiw__wptag(o, "IEND"); - stbiw__wpcrc(&o,0); - - STBIW_ASSERT(o == out + *out_len); - - return out; -} - -#ifndef STBI_WRITE_NO_STDIO -STBIWDEF int stbi_write_png(char const *filename, int x, int y, int comp, const void *data, int stride_bytes) -{ - FILE *f; - int len; - unsigned char *png = stbi_write_png_to_mem((unsigned char *) data, stride_bytes, x, y, comp, &len); - if (png == NULL) return 0; -#ifdef STBI_MSC_SECURE_CRT - if (fopen_s(&f, filename, "wb")) - f = NULL; -#else - f = fopen(filename, "wb"); -#endif - if (!f) { STBIW_FREE(png); return 0; } - fwrite(png, 1, len, f); - fclose(f); - STBIW_FREE(png); - return 1; -} -#endif - -STBIWDEF int stbi_write_png_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data, int stride_bytes) -{ - int len; - unsigned char *png = stbi_write_png_to_mem((unsigned char *) data, stride_bytes, x, y, comp, &len); - if (png == NULL) return 0; - func(context, png, len); - STBIW_FREE(png); - return 1; -} - - -/* *************************************************************************** - * - * JPEG writer - * - * This is based on Jon Olick's jo_jpeg.cpp: - * public domain Simple, Minimalistic JPEG writer - http://www.jonolick.com/code.html - */ - -static const unsigned char stbiw__jpg_ZigZag[] = { 0,1,5,6,14,15,27,28,2,4,7,13,16,26,29,42,3,8,12,17,25,30,41,43,9,11,18, - 24,31,40,44,53,10,19,23,32,39,45,52,54,20,22,33,38,46,51,55,60,21,34,37,47,50,56,59,61,35,36,48,49,57,58,62,63 }; - -static void stbiw__jpg_writeBits(stbi__write_context *s, int *bitBufP, int *bitCntP, const unsigned short *bs) { - int bitBuf = *bitBufP, bitCnt = *bitCntP; - bitCnt += bs[1]; - bitBuf |= bs[0] << (24 - bitCnt); - while(bitCnt >= 8) { - unsigned char c = (bitBuf >> 16) & 255; - stbiw__putc(s, c); - if(c == 255) { - stbiw__putc(s, 0); - } - bitBuf <<= 8; - bitCnt -= 8; - } - *bitBufP = bitBuf; - *bitCntP = bitCnt; -} - -static void stbiw__jpg_DCT(float *d0p, float *d1p, float *d2p, float *d3p, float *d4p, float *d5p, float *d6p, float *d7p) { - float d0 = *d0p, d1 = *d1p, d2 = *d2p, d3 = *d3p, d4 = *d4p, d5 = *d5p, d6 = *d6p, d7 = *d7p; - float z1, z2, z3, z4, z5, z11, z13; - - float tmp0 = d0 + d7; - float tmp7 = d0 - d7; - float tmp1 = d1 + d6; - float tmp6 = d1 - d6; - float tmp2 = d2 + d5; - float tmp5 = d2 - d5; - float tmp3 = d3 + d4; - float tmp4 = d3 - d4; - - // Even part - float tmp10 = tmp0 + tmp3; // phase 2 - float tmp13 = tmp0 - tmp3; - float tmp11 = tmp1 + tmp2; - float tmp12 = tmp1 - tmp2; - - d0 = tmp10 + tmp11; // phase 3 - d4 = tmp10 - tmp11; - - z1 = (tmp12 + tmp13) * 0.707106781f; // c4 - d2 = tmp13 + z1; // phase 5 - d6 = tmp13 - z1; - - // Odd part - tmp10 = tmp4 + tmp5; // phase 2 - tmp11 = tmp5 + tmp6; - tmp12 = tmp6 + tmp7; - - // The rotator is modified from fig 4-8 to avoid extra negations. - z5 = (tmp10 - tmp12) * 0.382683433f; // c6 - z2 = tmp10 * 0.541196100f + z5; // c2-c6 - z4 = tmp12 * 1.306562965f + z5; // c2+c6 - z3 = tmp11 * 0.707106781f; // c4 - - z11 = tmp7 + z3; // phase 5 - z13 = tmp7 - z3; - - *d5p = z13 + z2; // phase 6 - *d3p = z13 - z2; - *d1p = z11 + z4; - *d7p = z11 - z4; - - *d0p = d0; *d2p = d2; *d4p = d4; *d6p = d6; -} - -static void stbiw__jpg_calcBits(int val, unsigned short bits[2]) { - int tmp1 = val < 0 ? -val : val; - val = val < 0 ? val-1 : val; - bits[1] = 1; - while(tmp1 >>= 1) { - ++bits[1]; - } - bits[0] = val & ((1<0)&&(DU[end0pos]==0); --end0pos) { - } - // end0pos = first element in reverse order !=0 - if(end0pos == 0) { - stbiw__jpg_writeBits(s, bitBuf, bitCnt, EOB); - return DU[0]; - } - for(i = 1; i <= end0pos; ++i) { - int startpos = i; - int nrzeroes; - unsigned short bits[2]; - for (; DU[i]==0 && i<=end0pos; ++i) { - } - nrzeroes = i-startpos; - if ( nrzeroes >= 16 ) { - int lng = nrzeroes>>4; - int nrmarker; - for (nrmarker=1; nrmarker <= lng; ++nrmarker) - stbiw__jpg_writeBits(s, bitBuf, bitCnt, M16zeroes); - nrzeroes &= 15; - } - stbiw__jpg_calcBits(DU[i], bits); - stbiw__jpg_writeBits(s, bitBuf, bitCnt, HTAC[(nrzeroes<<4)+bits[1]]); - stbiw__jpg_writeBits(s, bitBuf, bitCnt, bits); - } - if(end0pos != 63) { - stbiw__jpg_writeBits(s, bitBuf, bitCnt, EOB); - } - return DU[0]; -} - -static int stbi_write_jpg_core(stbi__write_context *s, int width, int height, int comp, const void* data, int quality) { - // Constants that don't pollute global namespace - static const unsigned char std_dc_luminance_nrcodes[] = {0,0,1,5,1,1,1,1,1,1,0,0,0,0,0,0,0}; - static const unsigned char std_dc_luminance_values[] = {0,1,2,3,4,5,6,7,8,9,10,11}; - static const unsigned char std_ac_luminance_nrcodes[] = {0,0,2,1,3,3,2,4,3,5,5,4,4,0,0,1,0x7d}; - static const unsigned char std_ac_luminance_values[] = { - 0x01,0x02,0x03,0x00,0x04,0x11,0x05,0x12,0x21,0x31,0x41,0x06,0x13,0x51,0x61,0x07,0x22,0x71,0x14,0x32,0x81,0x91,0xa1,0x08, - 0x23,0x42,0xb1,0xc1,0x15,0x52,0xd1,0xf0,0x24,0x33,0x62,0x72,0x82,0x09,0x0a,0x16,0x17,0x18,0x19,0x1a,0x25,0x26,0x27,0x28, - 0x29,0x2a,0x34,0x35,0x36,0x37,0x38,0x39,0x3a,0x43,0x44,0x45,0x46,0x47,0x48,0x49,0x4a,0x53,0x54,0x55,0x56,0x57,0x58,0x59, - 0x5a,0x63,0x64,0x65,0x66,0x67,0x68,0x69,0x6a,0x73,0x74,0x75,0x76,0x77,0x78,0x79,0x7a,0x83,0x84,0x85,0x86,0x87,0x88,0x89, - 0x8a,0x92,0x93,0x94,0x95,0x96,0x97,0x98,0x99,0x9a,0xa2,0xa3,0xa4,0xa5,0xa6,0xa7,0xa8,0xa9,0xaa,0xb2,0xb3,0xb4,0xb5,0xb6, - 0xb7,0xb8,0xb9,0xba,0xc2,0xc3,0xc4,0xc5,0xc6,0xc7,0xc8,0xc9,0xca,0xd2,0xd3,0xd4,0xd5,0xd6,0xd7,0xd8,0xd9,0xda,0xe1,0xe2, - 0xe3,0xe4,0xe5,0xe6,0xe7,0xe8,0xe9,0xea,0xf1,0xf2,0xf3,0xf4,0xf5,0xf6,0xf7,0xf8,0xf9,0xfa - }; - static const unsigned char std_dc_chrominance_nrcodes[] = {0,0,3,1,1,1,1,1,1,1,1,1,0,0,0,0,0}; - static const unsigned char std_dc_chrominance_values[] = {0,1,2,3,4,5,6,7,8,9,10,11}; - static const unsigned char std_ac_chrominance_nrcodes[] = {0,0,2,1,2,4,4,3,4,7,5,4,4,0,1,2,0x77}; - static const unsigned char std_ac_chrominance_values[] = { - 0x00,0x01,0x02,0x03,0x11,0x04,0x05,0x21,0x31,0x06,0x12,0x41,0x51,0x07,0x61,0x71,0x13,0x22,0x32,0x81,0x08,0x14,0x42,0x91, - 0xa1,0xb1,0xc1,0x09,0x23,0x33,0x52,0xf0,0x15,0x62,0x72,0xd1,0x0a,0x16,0x24,0x34,0xe1,0x25,0xf1,0x17,0x18,0x19,0x1a,0x26, - 0x27,0x28,0x29,0x2a,0x35,0x36,0x37,0x38,0x39,0x3a,0x43,0x44,0x45,0x46,0x47,0x48,0x49,0x4a,0x53,0x54,0x55,0x56,0x57,0x58, - 0x59,0x5a,0x63,0x64,0x65,0x66,0x67,0x68,0x69,0x6a,0x73,0x74,0x75,0x76,0x77,0x78,0x79,0x7a,0x82,0x83,0x84,0x85,0x86,0x87, - 0x88,0x89,0x8a,0x92,0x93,0x94,0x95,0x96,0x97,0x98,0x99,0x9a,0xa2,0xa3,0xa4,0xa5,0xa6,0xa7,0xa8,0xa9,0xaa,0xb2,0xb3,0xb4, - 0xb5,0xb6,0xb7,0xb8,0xb9,0xba,0xc2,0xc3,0xc4,0xc5,0xc6,0xc7,0xc8,0xc9,0xca,0xd2,0xd3,0xd4,0xd5,0xd6,0xd7,0xd8,0xd9,0xda, - 0xe2,0xe3,0xe4,0xe5,0xe6,0xe7,0xe8,0xe9,0xea,0xf2,0xf3,0xf4,0xf5,0xf6,0xf7,0xf8,0xf9,0xfa - }; - // Huffman tables - static const unsigned short YDC_HT[256][2] = { {0,2},{2,3},{3,3},{4,3},{5,3},{6,3},{14,4},{30,5},{62,6},{126,7},{254,8},{510,9}}; - static const unsigned short UVDC_HT[256][2] = { {0,2},{1,2},{2,2},{6,3},{14,4},{30,5},{62,6},{126,7},{254,8},{510,9},{1022,10},{2046,11}}; - static const unsigned short YAC_HT[256][2] = { - {10,4},{0,2},{1,2},{4,3},{11,4},{26,5},{120,7},{248,8},{1014,10},{65410,16},{65411,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {12,4},{27,5},{121,7},{502,9},{2038,11},{65412,16},{65413,16},{65414,16},{65415,16},{65416,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {28,5},{249,8},{1015,10},{4084,12},{65417,16},{65418,16},{65419,16},{65420,16},{65421,16},{65422,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {58,6},{503,9},{4085,12},{65423,16},{65424,16},{65425,16},{65426,16},{65427,16},{65428,16},{65429,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {59,6},{1016,10},{65430,16},{65431,16},{65432,16},{65433,16},{65434,16},{65435,16},{65436,16},{65437,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {122,7},{2039,11},{65438,16},{65439,16},{65440,16},{65441,16},{65442,16},{65443,16},{65444,16},{65445,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {123,7},{4086,12},{65446,16},{65447,16},{65448,16},{65449,16},{65450,16},{65451,16},{65452,16},{65453,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {250,8},{4087,12},{65454,16},{65455,16},{65456,16},{65457,16},{65458,16},{65459,16},{65460,16},{65461,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {504,9},{32704,15},{65462,16},{65463,16},{65464,16},{65465,16},{65466,16},{65467,16},{65468,16},{65469,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {505,9},{65470,16},{65471,16},{65472,16},{65473,16},{65474,16},{65475,16},{65476,16},{65477,16},{65478,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {506,9},{65479,16},{65480,16},{65481,16},{65482,16},{65483,16},{65484,16},{65485,16},{65486,16},{65487,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {1017,10},{65488,16},{65489,16},{65490,16},{65491,16},{65492,16},{65493,16},{65494,16},{65495,16},{65496,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {1018,10},{65497,16},{65498,16},{65499,16},{65500,16},{65501,16},{65502,16},{65503,16},{65504,16},{65505,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {2040,11},{65506,16},{65507,16},{65508,16},{65509,16},{65510,16},{65511,16},{65512,16},{65513,16},{65514,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {65515,16},{65516,16},{65517,16},{65518,16},{65519,16},{65520,16},{65521,16},{65522,16},{65523,16},{65524,16},{0,0},{0,0},{0,0},{0,0},{0,0}, - {2041,11},{65525,16},{65526,16},{65527,16},{65528,16},{65529,16},{65530,16},{65531,16},{65532,16},{65533,16},{65534,16},{0,0},{0,0},{0,0},{0,0},{0,0} - }; - static const unsigned short UVAC_HT[256][2] = { - {0,2},{1,2},{4,3},{10,4},{24,5},{25,5},{56,6},{120,7},{500,9},{1014,10},{4084,12},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {11,4},{57,6},{246,8},{501,9},{2038,11},{4085,12},{65416,16},{65417,16},{65418,16},{65419,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {26,5},{247,8},{1015,10},{4086,12},{32706,15},{65420,16},{65421,16},{65422,16},{65423,16},{65424,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {27,5},{248,8},{1016,10},{4087,12},{65425,16},{65426,16},{65427,16},{65428,16},{65429,16},{65430,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {58,6},{502,9},{65431,16},{65432,16},{65433,16},{65434,16},{65435,16},{65436,16},{65437,16},{65438,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {59,6},{1017,10},{65439,16},{65440,16},{65441,16},{65442,16},{65443,16},{65444,16},{65445,16},{65446,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {121,7},{2039,11},{65447,16},{65448,16},{65449,16},{65450,16},{65451,16},{65452,16},{65453,16},{65454,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {122,7},{2040,11},{65455,16},{65456,16},{65457,16},{65458,16},{65459,16},{65460,16},{65461,16},{65462,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {249,8},{65463,16},{65464,16},{65465,16},{65466,16},{65467,16},{65468,16},{65469,16},{65470,16},{65471,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {503,9},{65472,16},{65473,16},{65474,16},{65475,16},{65476,16},{65477,16},{65478,16},{65479,16},{65480,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {504,9},{65481,16},{65482,16},{65483,16},{65484,16},{65485,16},{65486,16},{65487,16},{65488,16},{65489,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {505,9},{65490,16},{65491,16},{65492,16},{65493,16},{65494,16},{65495,16},{65496,16},{65497,16},{65498,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {506,9},{65499,16},{65500,16},{65501,16},{65502,16},{65503,16},{65504,16},{65505,16},{65506,16},{65507,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {2041,11},{65508,16},{65509,16},{65510,16},{65511,16},{65512,16},{65513,16},{65514,16},{65515,16},{65516,16},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}, - {16352,14},{65517,16},{65518,16},{65519,16},{65520,16},{65521,16},{65522,16},{65523,16},{65524,16},{65525,16},{0,0},{0,0},{0,0},{0,0},{0,0}, - {1018,10},{32707,15},{65526,16},{65527,16},{65528,16},{65529,16},{65530,16},{65531,16},{65532,16},{65533,16},{65534,16},{0,0},{0,0},{0,0},{0,0},{0,0} - }; - static const int YQT[] = {16,11,10,16,24,40,51,61,12,12,14,19,26,58,60,55,14,13,16,24,40,57,69,56,14,17,22,29,51,87,80,62,18,22, - 37,56,68,109,103,77,24,35,55,64,81,104,113,92,49,64,78,87,103,121,120,101,72,92,95,98,112,100,103,99}; - static const int UVQT[] = {17,18,24,47,99,99,99,99,18,21,26,66,99,99,99,99,24,26,56,99,99,99,99,99,47,66,99,99,99,99,99,99, - 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99}; - static const float aasf[] = { 1.0f * 2.828427125f, 1.387039845f * 2.828427125f, 1.306562965f * 2.828427125f, 1.175875602f * 2.828427125f, - 1.0f * 2.828427125f, 0.785694958f * 2.828427125f, 0.541196100f * 2.828427125f, 0.275899379f * 2.828427125f }; - - int row, col, i, k; - float fdtbl_Y[64], fdtbl_UV[64]; - unsigned char YTable[64], UVTable[64]; - - if(!data || !width || !height || comp > 4 || comp < 1) { - return 0; - } - - quality = quality ? quality : 90; - quality = quality < 1 ? 1 : quality > 100 ? 100 : quality; - quality = quality < 50 ? 5000 / quality : 200 - quality * 2; - - for(i = 0; i < 64; ++i) { - int uvti, yti = (YQT[i]*quality+50)/100; - YTable[stbiw__jpg_ZigZag[i]] = (unsigned char) (yti < 1 ? 1 : yti > 255 ? 255 : yti); - uvti = (UVQT[i]*quality+50)/100; - UVTable[stbiw__jpg_ZigZag[i]] = (unsigned char) (uvti < 1 ? 1 : uvti > 255 ? 255 : uvti); - } - - for(row = 0, k = 0; row < 8; ++row) { - for(col = 0; col < 8; ++col, ++k) { - fdtbl_Y[k] = 1 / (YTable [stbiw__jpg_ZigZag[k]] * aasf[row] * aasf[col]); - fdtbl_UV[k] = 1 / (UVTable[stbiw__jpg_ZigZag[k]] * aasf[row] * aasf[col]); - } - } - - // Write Headers - { - static const unsigned char head0[] = { 0xFF,0xD8,0xFF,0xE0,0,0x10,'J','F','I','F',0,1,1,0,0,1,0,1,0,0,0xFF,0xDB,0,0x84,0 }; - static const unsigned char head2[] = { 0xFF,0xDA,0,0xC,3,1,0,2,0x11,3,0x11,0,0x3F,0 }; - const unsigned char head1[] = { 0xFF,0xC0,0,0x11,8,(unsigned char)(height>>8),STBIW_UCHAR(height),(unsigned char)(width>>8),STBIW_UCHAR(width), - 3,1,0x11,0,2,0x11,1,3,0x11,1,0xFF,0xC4,0x01,0xA2,0 }; - s->func(s->context, (void*)head0, sizeof(head0)); - s->func(s->context, (void*)YTable, sizeof(YTable)); - stbiw__putc(s, 1); - s->func(s->context, UVTable, sizeof(UVTable)); - s->func(s->context, (void*)head1, sizeof(head1)); - s->func(s->context, (void*)(std_dc_luminance_nrcodes+1), sizeof(std_dc_luminance_nrcodes)-1); - s->func(s->context, (void*)std_dc_luminance_values, sizeof(std_dc_luminance_values)); - stbiw__putc(s, 0x10); // HTYACinfo - s->func(s->context, (void*)(std_ac_luminance_nrcodes+1), sizeof(std_ac_luminance_nrcodes)-1); - s->func(s->context, (void*)std_ac_luminance_values, sizeof(std_ac_luminance_values)); - stbiw__putc(s, 1); // HTUDCinfo - s->func(s->context, (void*)(std_dc_chrominance_nrcodes+1), sizeof(std_dc_chrominance_nrcodes)-1); - s->func(s->context, (void*)std_dc_chrominance_values, sizeof(std_dc_chrominance_values)); - stbiw__putc(s, 0x11); // HTUACinfo - s->func(s->context, (void*)(std_ac_chrominance_nrcodes+1), sizeof(std_ac_chrominance_nrcodes)-1); - s->func(s->context, (void*)std_ac_chrominance_values, sizeof(std_ac_chrominance_values)); - s->func(s->context, (void*)head2, sizeof(head2)); - } - - // Encode 8x8 macroblocks - { - static const unsigned short fillBits[] = {0x7F, 7}; - const unsigned char *imageData = (const unsigned char *)data; - int DCY=0, DCU=0, DCV=0; - int bitBuf=0, bitCnt=0; - // comp == 2 is grey+alpha (alpha is ignored) - int ofsG = comp > 2 ? 1 : 0, ofsB = comp > 2 ? 2 : 0; - int x, y, pos; - for(y = 0; y < height; y += 8) { - for(x = 0; x < width; x += 8) { - float YDU[64], UDU[64], VDU[64]; - for(row = y, pos = 0; row < y+8; ++row) { - for(col = x; col < x+8; ++col, ++pos) { - int p = (stbi__flip_vertically_on_write ? height-1-row : row)*width*comp + col*comp; - float r, g, b; - if(row >= height) { - p -= width*comp*(row+1 - height); - } - if(col >= width) { - p -= comp*(col+1 - width); - } - - r = imageData[p+0]; - g = imageData[p+ofsG]; - b = imageData[p+ofsB]; - YDU[pos]=+0.29900f*r+0.58700f*g+0.11400f*b-128; - UDU[pos]=-0.16874f*r-0.33126f*g+0.50000f*b; - VDU[pos]=+0.50000f*r-0.41869f*g-0.08131f*b; - } - } - - DCY = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, YDU, fdtbl_Y, DCY, YDC_HT, YAC_HT); - DCU = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, UDU, fdtbl_UV, DCU, UVDC_HT, UVAC_HT); - DCV = stbiw__jpg_processDU(s, &bitBuf, &bitCnt, VDU, fdtbl_UV, DCV, UVDC_HT, UVAC_HT); - } - } - - // Do the bit alignment of the EOI marker - stbiw__jpg_writeBits(s, &bitBuf, &bitCnt, fillBits); - } - - // EOI - stbiw__putc(s, 0xFF); - stbiw__putc(s, 0xD9); - - return 1; -} - -STBIWDEF int stbi_write_jpg_to_func(stbi_write_func *func, void *context, int x, int y, int comp, const void *data, int quality) -{ - stbi__write_context s; - stbi__start_write_callbacks(&s, func, context); - return stbi_write_jpg_core(&s, x, y, comp, (void *) data, quality); -} - - -#ifndef STBI_WRITE_NO_STDIO -STBIWDEF int stbi_write_jpg(char const *filename, int x, int y, int comp, const void *data, int quality) -{ - stbi__write_context s; - if (stbi__start_write_file(&s,filename)) { - int r = stbi_write_jpg_core(&s, x, y, comp, data, quality); - stbi__end_write_file(&s); - return r; - } else - return 0; -} -#endif - -#endif // STB_IMAGE_WRITE_IMPLEMENTATION - -/* Revision history - 1.09 (2018-02-11) - fix typo in zlib quality API, improve STB_I_W_STATIC in C++ - 1.08 (2018-01-29) - add stbi__flip_vertically_on_write, external zlib, zlib quality, choose PNG filter - 1.07 (2017-07-24) - doc fix - 1.06 (2017-07-23) - writing JPEG (using Jon Olick's code) - 1.05 ??? - 1.04 (2017-03-03) - monochrome BMP expansion - 1.03 ??? - 1.02 (2016-04-02) - avoid allocating large structures on the stack - 1.01 (2016-01-16) - STBIW_REALLOC_SIZED: support allocators with no realloc support - avoid race-condition in crc initialization - minor compile issues - 1.00 (2015-09-14) - installable file IO function - 0.99 (2015-09-13) - warning fixes; TGA rle support - 0.98 (2015-04-08) - added STBIW_MALLOC, STBIW_ASSERT etc - 0.97 (2015-01-18) - fixed HDR asserts, rewrote HDR rle logic - 0.96 (2015-01-17) - add HDR output - fix monochrome BMP - 0.95 (2014-08-17) - add monochrome TGA output - 0.94 (2014-05-31) - rename private functions to avoid conflicts with stb_image.h - 0.93 (2014-05-27) - warning fixes - 0.92 (2010-08-01) - casts to unsigned char to fix warnings - 0.91 (2010-07-17) - first public release - 0.90 first internal release -*/ - -/* ------------------------------------------------------------------------------- -This software is available under 2 licenses -- choose whichever you prefer. ------------------------------------------------------------------------------- -ALTERNATIVE A - MIT License -Copyright (c) 2017 Sean Barrett -Permission is hereby granted, free of charge, to any person obtaining a copy of -this software and associated documentation files (the "Software"), to deal in -the Software without restriction, including without limitation the rights to -use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies -of the Software, and to permit persons to whom the Software is furnished to do -so, subject to the following conditions: -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. ------------------------------------------------------------------------------- -ALTERNATIVE B - Public Domain (www.unlicense.org) -This is free and unencumbered software released into the public domain. -Anyone is free to copy, modify, publish, use, compile, sell, or distribute this -software, either in source code form or as a compiled binary, for any purpose, -commercial or non-commercial, and by any means. -In jurisdictions that recognize copyright laws, the author or authors of this -software dedicate any and all copyright interest in the software to the public -domain. We make this dedication for the benefit of the public at large and to -the detriment of our heirs and successors. We intend this dedication to be an -overt act of relinquishment in perpetuity of all present and future rights to -this software under copyright law. -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN -ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION -WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ------------------------------------------------------------------------------- -*/ diff --git a/util.h b/util.h deleted file mode 100644 index 821c3e6..0000000 --- a/util.h +++ /dev/null @@ -1,49 +0,0 @@ -#ifndef UTIL_H -#define UTIL_H - -#include // fprintf() -#include // exit() - -#define MIN(a, b) ((a) < (b) ? (a) : (b)) -#define MAX(a, b) ((a) > (b) ? (a) : (b)) -#define UNUSED(a) ((void) (a)) - -#define D(...) do { \ - fprintf(stderr, "DEBUG: %s(): ", __func__); \ - fprintf(stderr, __VA_ARGS__); \ - fputs("\n", stderr); \ -} while (0) - -#define die(...) do { \ - fputs("FATAL: ", stderr); \ - fprintf(stderr, __VA_ARGS__); \ - fputs("\n", stderr); \ - exit(EXIT_FAILURE); \ -} while (0) - -// calculate squared euclidean distance between two points -static inline float -distance_squared( - const size_t num_floats, - const float * const a, - const float * const b -) { - float r = 0.0; - - for (size_t i = 0; i < num_floats; i++) { - r += (b[i] - a[i]) * (b[i] - a[i]); - } - - // return squared distance - return r; -} - -static inline float -km_score( - const float mean_distance, - const size_t num_empty -) { - return 1.0 / (mean_distance + num_empty); -} - -#endif /* UTIL_H */ -- cgit v1.2.3