#include <gdk-pixbuf/gdk-pixbuf.h>
#include <glib.h>
#include <string.h>

#ifndef INTERVAL
#define INTERVAL 8
#endif
#define BLOCKS (256 / INTERVAL)

struct cell {
	int count;
	struct cell *max;
};

struct hst {
	struct cell cells[BLOCKS][BLOCKS][BLOCKS];
	GList *list;
};

static inline void
pixbuf_get(GdkPixbuf *pixbuf, guchar **pixels,
	   gint *width, gint *height, gint *nch, gint *stride)
{
	if (pixels) *pixels = gdk_pixbuf_get_pixels(pixbuf);
	if (width)  *width  = gdk_pixbuf_get_width(pixbuf);
	if (height) *height = gdk_pixbuf_get_height(pixbuf);
	if (nch)    *nch    = gdk_pixbuf_get_n_channels(pixbuf);
	if (stride) *stride = gdk_pixbuf_get_rowstride(pixbuf);
}

static gint
compare_cell(gconstpointer a, gconstpointer b)
{
	const struct cell *cell1 = a;
	const struct cell *cell2 = b;

	return cell1->count - cell2->count;
}

static void
calculate_index(struct hst *hst, struct cell *cell, int *i, int *j, int *k)
{
	size_t offset = cell - &hst->cells[0][0][0];

	*i = offset / (BLOCKS*BLOCKS);
	*j = (offset - *i * (BLOCKS*BLOCKS)) / BLOCKS;
	*k = (offset - *i * (BLOCKS*BLOCKS) - *j * BLOCKS);
}

static void
add_neighbors(struct hst *hst, struct cell *cell)
{
	int i, j, k, ci, cj, ck;

	calculate_index(hst, cell, &ci, &cj, &ck);
	for (i = MAX(ci - 1, 0); i < MIN(ci + 1, BLOCKS); ++i) {
		for (j = MAX(cj - 1, 0); j < MIN(cj + 1, BLOCKS); ++j) {
			for (k = MAX(ck - 1, 0); k < MIN(ck + 1, BLOCKS); ++k) {
				if (hst->cells[i][j][k].count <= cell->count &&
				    hst->cells[i][j][k].max == NULL) {
				    	hst->cells[i][j][k].max = cell->max;
					add_neighbors(hst, &hst->cells[i][j][k]);
				}
			}
		}
	}
	hst->list = g_list_remove(hst->list, cell);
}

int
main(int argc, char *argv[])
{
	GdkPixbuf *image;
	gint x, y, width, height, nch, stride;
	int i, j, k;
	guchar *p, *row;
	int cluster = 0;
	struct hst hst;
	struct cell *cell;

	g_type_init();
	memset(&hst, 0, sizeof hst);

	if (argc < 2)
		return 1;

	image = gdk_pixbuf_new_from_file(argv[1], NULL);
	if (!image)
		return 1;

	pixbuf_get(image, &row, &width, &height, &nch, &stride);
	for (y = 0; y < height; ++y, row += stride)
		for (x = 0, p = row; x < width; ++x, p += nch)
			++hst.cells[p[0] / INTERVAL][p[1] / INTERVAL][p[2] / INTERVAL].count;

	for (i = 0; i < BLOCKS; ++i)
		for (j = 0; j < BLOCKS; ++j)
			for (k = 0; k < BLOCKS; ++k)
				hst.list = g_list_insert_sorted(hst.list,
								&hst.cells[i][j][k],
								compare_cell);

	while (hst.list) {
		cell = g_list_nth_data(hst.list, 0);
		cell->max = cell;
		add_neighbors(&hst, cell);
		++cluster;
	}
	printf("cluster: %d\n", cluster);

	pixbuf_get(image, &row, &width, &height, &nch, &stride);
	for (y = 0; y < height; ++y, row += stride) {
		for (x = 0, p = row; x < width; ++x, p += nch) {
			cell = &hst.cells[p[0] / INTERVAL][p[1] / INTERVAL][p[2] / INTERVAL];

			calculate_index(&hst, cell->max, &i, &j, &k);
			p[0] = i * INTERVAL + INTERVAL / 2;
			p[1] = j * INTERVAL + INTERVAL / 2;
			p[2] = k * INTERVAL + INTERVAL / 2;
		}
	}

	gdk_pixbuf_save(image, argc > 2 ? argv[2] : "hst.tiff", "tiff", NULL, NULL);

	return 0;
}