aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/km-draw.c31
-rw-r--r--src/km-find.c115
-rw-r--r--src/km-init-forgy.c51
-rw-r--r--src/km-init-kmeans.c124
-rw-r--r--src/km-init-rand.c43
-rw-r--r--src/km-init.c85
-rw-r--r--src/km-load.c112
-rw-r--r--src/km-rand-erand48.c89
-rw-r--r--src/km-rand-libc.c56
-rw-r--r--src/km-rand-path.c152
-rw-r--r--src/km-rand.c33
-rw-r--r--src/km-set.c267
-rw-r--r--src/km-solve.c200
-rw-r--r--src/km.h230
-rw-r--r--src/main.c398
-rw-r--r--src/stb_image_write.h1568
-rw-r--r--src/util.h49
17 files changed, 3603 insertions, 0 deletions
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 <stdint.h> // 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 <stdbool.h> // bool
+#include <math.h> // 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 <stdbool.h> // bool
+#include <string.h> // 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 <stdbool.h> // bool
+#include <string.h> // memset()
+#include <float.h> // 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 <stdbool.h> // bool
+#include <string.h> // 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 <stdbool.h> // bool
+#include <stddef.h> // size_t
+#include <string.h> // 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 <stdbool.h> // bool
+#include <stdio.h> // fscanf()
+#include <string.h> // strerror()
+#include <errno.h> // 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 <stdbool.h> // bool
+#define _DEFAULT_SOURCE
+#include <stdlib.h> // drand48_r()
+#include <stdio.h> // 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 <stdbool.h> // 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 <stdbool.h> // bool
+#include <string.h> // memset()
+#include <stdio.h> // 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 <stdbool.h>
+#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 <stdbool.h> // bool
+#include <stdint.h> // size_t
+#include <string.h> // memcpy()
+#include <stdlib.h> // rand()
+#include <errno.h> // 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 <stdbool.h> // bool
+#include <string.h> // memset()
+#include <float.h> // FLT_MAX
+#include <math.h> // 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 <stddef.h> // size_t
+#include <stdint.h> // uint8_t
+#include <stdio.h> // 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 <stdbool.h> // bool
+#include <stdio.h> // fprintf()
+#include <unistd.h> // EXIT_{FAILURE,SUCCESS}
+#include <stdlib.h> // exit()
+#include <string.h> // memset()
+#include <math.h> // 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] <png_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 <stdio.h>
+#endif // STBI_WRITE_NO_STDIO
+
+#include <stdarg.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#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 <assert.h>
+#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<<bits[1])-1);
+}
+
+static int stbiw__jpg_processDU(stbi__write_context *s, int *bitBuf, int *bitCnt, float *CDU, float *fdtbl, int DC, const unsigned short HTDC[256][2], const unsigned short HTAC[256][2]) {
+ const unsigned short EOB[2] = { HTAC[0x00][0], HTAC[0x00][1] };
+ const unsigned short M16zeroes[2] = { HTAC[0xF0][0], HTAC[0xF0][1] };
+ int dataOff, i, diff, end0pos;
+ int DU[64];
+
+ // DCT rows
+ for(dataOff=0; dataOff<64; dataOff+=8) {
+ stbiw__jpg_DCT(&CDU[dataOff], &CDU[dataOff+1], &CDU[dataOff+2], &CDU[dataOff+3], &CDU[dataOff+4], &CDU[dataOff+5], &CDU[dataOff+6], &CDU[dataOff+7]);
+ }
+ // DCT columns
+ for(dataOff=0; dataOff<8; ++dataOff) {
+ stbiw__jpg_DCT(&CDU[dataOff], &CDU[dataOff+8], &CDU[dataOff+16], &CDU[dataOff+24], &CDU[dataOff+32], &CDU[dataOff+40], &CDU[dataOff+48], &CDU[dataOff+56]);
+ }
+ // Quantize/descale/zigzag the coefficients
+ for(i=0; i<64; ++i) {
+ float v = CDU[i]*fdtbl[i];
+ // DU[stbiw__jpg_ZigZag[i]] = (int)(v < 0 ? ceilf(v - 0.5f) : floorf(v + 0.5f));
+ // ceilf() and floorf() are C99, not C89, but I /think/ they're not needed here anyway?
+ DU[stbiw__jpg_ZigZag[i]] = (int)(v < 0 ? v - 0.5f : v + 0.5f);
+ }
+
+ // Encode DC
+ diff = DU[0] - DC;
+ if (diff == 0) {
+ stbiw__jpg_writeBits(s, bitBuf, bitCnt, HTDC[0]);
+ } else {
+ unsigned short bits[2];
+ stbiw__jpg_calcBits(diff, bits);
+ stbiw__jpg_writeBits(s, bitBuf, bitCnt, HTDC[bits[1]]);
+ stbiw__jpg_writeBits(s, bitBuf, bitCnt, bits);
+ }
+ // Encode ACs
+ end0pos = 63;
+ for(; (end0pos>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 <stdio.h> // fprintf()
+#include <stdlib.h> // 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 */
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-22 03:00:08 +0000