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

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

struct cell {
	int i, j, k;
	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
add_neighbors(struct hst *hst, struct cell *cell)
{
	int i, j, k;	

	for (i = MAX(cell->i - 1, 0); i < MIN(cell->i + 1, BLOCKS); ++i) {
		for (j = MAX(cell->j - 1, 0); j < MIN(cell->j + 1, BLOCKS); ++j) {
			for (k = MAX(cell->k - 1, 0); k < MIN(cell->k + 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();

	if (argc < 2)
		return 1;

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

	hst.list = NULL;
	for (i = 0; i < BLOCKS; ++i) {
		for (j = 0; j < BLOCKS; ++j) {
			for (k = 0; k < BLOCKS; ++k) {
				hst.cells[i][j][k].i = i;
				hst.cells[i][j][k].j = j;
				hst.cells[i][j][k].k = k;
				hst.cells[i][j][k].count = 0;
				hst.cells[i][j][k].max = NULL;
			}
		}
	}

	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];

			p[0] = cell->max->i * INTERVAL;
			p[1] = cell->max->j * INTERVAL;
			p[2] = cell->max->k * INTERVAL;
		}
	}

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

	return 0;
}