Add harris detector

2 files changed, 235 insertions, 1 deletions

diff --git a/Makefile b/Makefile
index 0e7e3a5..45ce1bf 100644
--- a/Makefile
+++ b/Makefile
@@ -25,7 +25,10 @@ cvg_LIBS=$(shell pkg-config --libs gtk+-3.0 glib-2.0 gmodule-2.0 opencv)
 hst_CFLAGS=$(shell pkg-config --cflags gdk-pixbuf-2.0 glib-2.0)
 hst_LIBS=$(shell pkg-config --libs gdk-pixbuf-2.0 glib-2.0)
 
-PROGS = wimmel wimmel_gl roi capture unvignette cvg hst
+harris_CFLAGS=$(shell pkg-config --cflags gdk-pixbuf-2.0 glib-2.0)
+harris_LIBS=$(shell pkg-config --libs gdk-pixbuf-2.0 glib-2.0)
+
+PROGS = wimmel wimmel_gl roi capture unvignette cvg hst harris
 OBJS = $(PROGS:=.o) util.o
 
 all: $(PROGS)
@@ -37,6 +40,7 @@ capture: capture.o
 unvignette: unvignette.o
 cvg: cvg.o
 hst: hst.o
+harris: harris.o
 
 .PHONY: clean all
 
diff --git a/harris.c b/harris.c
new file mode 100644
index 0000000..fab5141
--- /dev/null
+++ b/harris.c
@@ -0,0 +1,230 @@
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <gdk-pixbuf/gdk-pixbuf.h>
+
+struct mat2 {
+	float a1, a2, b1, b2;
+	float mc;
+};
+
+static inline float
+det(struct mat2 *m)
+{
+	return m->a1 * m->b2 - m->a2 * m->b1;
+}
+
+static inline float
+trace(struct mat2 *m)
+{
+	return m->a1 + m->b2;
+}
+
+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 void
+grey_scale(GdkPixbuf *image)
+{
+	gint Y, x, y, width, height, nch, stride;
+	guchar *p, *row;
+
+	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) {
+			Y = 0.299 * p[0] + 0.587 * p[1] + 0.115 * p[2];
+			p[0] = p[1] = p[2] = Y;
+		}
+	}
+}
+
+int *
+conv(GdkPixbuf *pixbuf, gint *operator, gint block_size)
+{
+	gint Y, i, j, x, y, width, height, nch, stride;
+	guchar *p, *row, *v;
+	int *out;
+	
+	pixbuf_get(pixbuf, &row, &width, &height, &nch, &stride);
+
+	out = malloc(width * height * sizeof *out);
+	if (!out)
+		return NULL;
+	for (y = 0; y < height; ++y, row += stride) {
+		for (x = 0, p = row; x < width; ++x, p += nch) {
+			Y = 0;
+			for (i = 0; i < block_size; ++i) {
+				for (j = 0; j < block_size; ++j) {
+					/* TODO: fix clamping */
+					v = p + (CLAMP((j-block_size/2), 0, width-1)  * stride +
+						 CLAMP((i-block_size/2), 0, height-1) * nch);
+					Y += v[0] * operator[j*block_size + i];
+				}
+			}
+			out[y * width + x] = Y/4;
+		}
+	}
+
+	return out;
+}
+
+struct mat2 *
+multiply(int *deriv_x, int *deriv_y, int width, int height)
+{
+	struct mat2 *m, *pm;
+	int *px, *py;
+	int i, n;
+
+	n = width * height;
+	m = calloc(n, sizeof *m);
+	if (!m)
+		return NULL;
+	px = deriv_x;
+	py = deriv_y;
+	pm = m;
+	for (i = 0; i < n; ++i) {
+		pm->a1          = *px * *px;
+		pm->a2 = pm->b1 = *px * *py;
+		pm->b2          = *py * *py;
+
+		++px; ++py; ++pm;
+	}
+
+	return m;
+}
+
+struct mat2 *
+calc_window_sum(struct mat2 *m, int width, int height, int wsize)
+{
+	struct mat2 *avg, *v, *r;
+	float w;
+	int x, y, i, j;
+
+	avg = calloc(width * height, sizeof *avg);
+	if (!avg)
+		return NULL;
+
+	r = avg;
+	for (y = 0; y < height; ++y) {
+		for (x = 0; x < width; ++x) {
+			for (i = -wsize/2; i < wsize/2; ++i) {
+				for (j = -wsize/2; j < wsize/2; ++j) {
+					v = &m[(CLAMP(y+i, 0, height-1) * width
+						+ CLAMP(x+j, 0, width-1))];
+					w = (1.0 / (2 * 3.14 * 0.16) *
+					     exp(- i*i - j*j / 0.16));
+
+					r->a1 += w * v->a1;
+					r->a2 += w * v->a2;
+					r->b1 += w * v->b1;
+					r->b2 += w * v->b2;
+				}
+			}
+			++r;
+		}
+	}
+
+	return avg;
+}
+
+gboolean
+check_neighbours(struct mat2 *m, int x, int y, int width, int height)
+{
+	int i, j;
+	gboolean max = TRUE;
+
+	for (j = MAX(0, y-1); j < MIN(y+1, height-1); ++j) {
+		for (i = MAX(0, x-1); i < MIN(x+1, width-1); ++i) {
+			if (m[j * width + i].mc > m[y * width + x].mc)
+				max = FALSE;
+		}
+	}
+
+	return max;
+}
+
+int
+main(int argc, char *argv[])
+{
+	GdkPixbuf *image, *grey;
+	int *deriv_x, *deriv_y;
+	gint sobel_x[] = {
+		1, 0, -1,
+		2, 0, -2,
+		1, 0, -1
+	};
+	gint sobel_y[] = {
+		1,  2,  1,
+		0,  0,  0,
+		-1, -2, -1
+	};
+	const float k = 0.04;
+	int x, y, width, height, size, nch, stride;
+	struct mat2 *m, *n, *c;
+	guchar *p, *row;
+	float max_mc, threshold;
+
+	if (argc < 2)
+		return 1;
+
+	g_type_init();
+	image = gdk_pixbuf_new_from_file(argv[1], NULL);
+	if (!image)
+		return 1;
+	width = gdk_pixbuf_get_width(image);
+	height = gdk_pixbuf_get_height(image);
+
+	grey = gdk_pixbuf_copy(image);
+	grey_scale(grey);
+
+	deriv_x = conv(grey, sobel_x, 3);
+	deriv_y = conv(grey, sobel_y, 3);
+	gdk_pixbuf_save(grey, argc > 3 ? argv[3] : "grey.tiff", "tiff",
+			NULL, NULL);
+
+	m = multiply(deriv_x, deriv_y, width, height);
+	n = calc_window_sum(m, width, height, 9);
+	free(m);
+
+	max_mc = 0;
+	size = width * height;
+	for (c = n; c < n+size; ++c) {
+		c->mc = det(c) - k * trace(c) * trace(c);
+		if (c->mc > max_mc)
+			max_mc = c->mc;
+	}
+	printf("max_mc: %f\n", max_mc);
+
+	threshold = 0.0001 * max_mc;
+	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) {
+			if (n[y * width + x].mc > threshold &&
+			    check_neighbours(n, x, y, width, height)) {
+				p[0] = 255;
+				p[1] = 0;
+				p[2] = 0;
+				printf("corner at: %d, %d\n", x, y);
+			}
+		}
+	}
+
+	gdk_pixbuf_save(image, argc > 2 ? argv[2] : "harris.tiff", "tiff",
+			NULL, NULL);
+
+	free(n);
+	free(deriv_x);
+	free(deriv_y);
+	g_object_unref(image);
+	g_object_unref(grey);
+
+	return 0;
+}