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

int
main(int argc, char *argv[])
{
	char *file;
	GdkPixbuf *input;
	GdkPixbuf *match_tmp, *match;
	int x, y, i, j;
	int x1, y1, mwidth, mheight, width, height;
	gint stride, nch;
	guchar *pix, *row;

	g_type_init();

	if (argc < 5)
		exit(EXIT_FAILURE);
	
	file = argv[1];
	x1 = atoi(argv[2]);
	y1 = atoi(argv[3]);
	mwidth = atoi(argv[4]);
	mheight = atoi(argv[5]);

	input = gdk_pixbuf_new_from_file(file, NULL);
	if (!input)
		exit(EXIT_FAILURE);

	match_tmp = gdk_pixbuf_new_subpixbuf(input, x1, y1, mwidth, mheight);
	match = gdk_pixbuf_copy(match_tmp);
	g_object_unref(match_tmp);

	struct match_interval {
		gint min, max;
	} interval[3];

	nch = gdk_pixbuf_get_n_channels(match);
	stride = gdk_pixbuf_get_rowstride(match);
	gint s = mwidth * mheight;

	gint middle[3] = { 0, 0, 0 };
	for (y = 0, row = gdk_pixbuf_get_pixels(match);
	     y < mheight; ++y, row += stride) {
		for (x = 0, pix = row; x < mwidth; ++x, pix += nch) {
			for (i = 0; i < 3; ++i)
				middle[i] += pix[i];
		}
	}
	for (i = 0; i < 3; ++i)
		middle[i] /= s;

	gint sigma[3] = { 0, 0, 0 };
	for (y = 0, row = gdk_pixbuf_get_pixels(match);
	     y < mheight; ++y, row += stride) {
		for (x = 0, pix = row; x < mwidth; ++x, pix += nch) {
			for (i = 0; i < 3; ++i)
				sigma[i] +=
					(pix[i] - middle[i]) *
					(pix[i] - middle[i]);
		}
	}

	for (i = 0; i < 3; ++i)
		sigma[i] = sqrt(sigma[i] / (s - 1));

	printf("middle: %3d, %3d, %3d\n", middle[0], middle[1], middle[2]);
	printf("sigma:  %3d, %3d, %3d\n", sigma[0], sigma[1], sigma[2]);

	for (i = 0; i < 3; ++i) {
		interval[i].min = middle[i] - 2 * sigma[i];
		interval[i].max = middle[i] + 2 * sigma[i];
		printf("interval %d: [%3d - %3d]\n",
		       i, interval[i].min, interval[i].max);
	}

	nch = gdk_pixbuf_get_n_channels(input);
	stride = gdk_pixbuf_get_rowstride(input);
	width = gdk_pixbuf_get_width(input);
	height = gdk_pixbuf_get_height(input);
	for (y = 0, row = gdk_pixbuf_get_pixels(input);
	     y < height; ++y, row += stride) {
		for (x = 0, pix = row; x < width; ++x, pix += nch) {
			int H = 1;
			for (i = 0; i < 3; ++i) {
				if (!(pix[i] >= interval[i].min &&
				      pix[i] <= interval[i].max))
				    H = 0;
			}
			if (H)
				pix[0] = pix[1] = pix[2] = 0;
		}
	}

	gdk_pixbuf_save(input, "output_roi.tiff", "tiff", NULL, NULL);
	g_object_unref(input);
	g_object_unref(match);
}