Add void cluster

7 files changed, 407 insertions, 22 deletions

diff --git a/Makefile b/Makefile
index ce169c7..0ba2746 100644
--- a/Makefile
+++ b/Makefile
@@ -1,8 +1,8 @@
-CFLAGS=-Wall
-LDFLAGS=-Wall -lm 
+CFLAGS=-Wall -g
+LDFLAGS=-Wall -lm -g
 GTKFLAGS=`pkg-config --cflags --libs gtk+-2.0`
 
-all: filters fft-test bluenoise trans noise2
+all: filters fft-test bluenoise trans noise2 voidcluster
 
 filters: fft.c fft.h filters.c
 	$(CC) -o filters fft.c filters.c $(LDFLAGS)
@@ -18,3 +18,6 @@ bluenoise: fft.c fft.h bluenoise.c gtk-utils.c
 
 noise2: noise2.c fft.c fft.h image.c image.h
 	$(CC) $(GTKFLAGS) -o noise2 noise2.c image.c fft.c $(LDFLAGS)
+
+voidcluster: voidcluster.c image.h fft.c
+	$(CC) $(GTKFLAGS) -o voidcluster fft.c voidcluster.c image.c $(LDFLAGS)
\ No newline at end of file
diff --git a/complex.h b/complex.h
index 58260a5..884656b 100644
--- a/complex.h
+++ b/complex.h
@@ -21,6 +21,15 @@ complex_mul (complex_t a, complex_t b)
 }
 
 static inline complex_t
+complex_smul (complex_t a, double s)
+{
+    a.re *= s;
+    a.im *= s;
+
+    return a;
+}
+
+static inline complex_t
 complex_add (complex_t a, complex_t b)
 {
     complex_t r;
diff --git a/image.c b/image.c
index 7b64513..a4ae29b 100644
--- a/image.c
+++ b/image.c
@@ -1,6 +1,7 @@
 #include <glib.h>
 #include <gtk/gtk.h>
 #include "image.h"
+#include "fft.h"
 
 complex_image_t *
 complex_image_new (int width, int height)
@@ -50,6 +51,27 @@ complex_image_subtract (complex_image_t *image, complex_image_t *other)
     }
 }
 
+void
+complex_image_mul (complex_image_t *image, double m)
+{
+    int i, j;
+    int h = image->height;
+    int w = image->width;
+    
+    for (i = 0; i < h; ++i)
+    {
+	for (j = 0; j < w; ++j)
+	{
+	    int idx = i * w + j;
+	    
+	    image->red[idx] = complex_smul (image->red[idx], m);
+	    image->green[idx] = complex_smul (image->green[idx], m);
+	    image->blue[idx] = complex_smul (image->blue[idx], m);
+	}
+    }
+}
+
+
 complex_image_t *
 complex_image_from_pixbuf (GdkPixbuf *pixbuf)
 {
@@ -215,3 +237,20 @@ complex_image_show (const char *title,
 
     gtk_main ();
 }
+
+void
+complex_image_fft (complex_image_t *image)
+{
+    fft_shift_2d (image->red, image->width);
+    fft_shift_2d (image->green, image->width);
+    fft_shift_2d (image->blue, image->width);
+}
+
+void
+complex_image_ifft (complex_image_t *image)
+{
+    ifft_shift_2d (image->red, image->width);
+    ifft_shift_2d (image->green, image->width);
+    ifft_shift_2d (image->blue, image->width);
+}
+
diff --git a/image.h b/image.h
index 1363eed..404e2f1 100644
--- a/image.h
+++ b/image.h
@@ -22,6 +22,9 @@ void
 complex_image_subtract (complex_image_t *image,
 			complex_image_t *other);
 
+void
+complex_image_mul (complex_image_t *image, double m);
+
 complex_image_t *
 complex_image_from_pixbuf (GdkPixbuf *pixbuf);
 
@@ -39,3 +42,9 @@ void
 complex_image_show (const char *title,
 		    complex_image_t *image,
 		    convert_type_t convert);
+void
+complex_image_fft (complex_image_t *image);
+
+void
+complex_image_ifft (complex_image_t *image);
+
diff --git a/noise2.c b/noise2.c
index b308ec2..606bd23 100644
--- a/noise2.c
+++ b/noise2.c
@@ -5,6 +5,7 @@
 #include <gdk/gdkpixbuf.h>
 #include "fft.h"
 #include "image.h"
+
 static void
 image_fft (complex_image_t *image)
 {
diff --git a/pngtrans.c b/pngtrans.c
index c1adf25..32df81f 100644
--- a/pngtrans.c
+++ b/pngtrans.c
@@ -6,29 +6,13 @@
 #include "fft.h"
 #include "image.h"
 
-static void
-image_fft (complex_image_t *image)
-{
-    fft_shift_2d (image->red, image->width);
-    fft_shift_2d (image->green, image->width);
-    fft_shift_2d (image->blue, image->width);
-}
-
-static void
-image_ifft (complex_image_t *image)
-{
-    ifft_shift_2d (image->red, image->width);
-    ifft_shift_2d (image->green, image->width);
-    ifft_shift_2d (image->blue, image->width);
-}
-
 static complex_t
 filter (complex_t c, double dist)
 {
     double m = complex_mag (c);
     double arg = complex_arg (c);
     
-    return complex_from_mag_arg (m * (1/(pow (sqrt (dist) / 16.0, 4) + 1)) , arg);
+    return complex_from_mag_arg (m * (1/(pow (sqrt (dist) / 128.0, 4) + 1)) , arg);
 }
 
 static void
@@ -82,7 +66,7 @@ main (int argc, char **argv)
     
     g_print ("width, height %d %d\n", image->width, image->height);
     
-    image_fft (image);
+    complex_image_fft (image);
     
     complex_image_show ("test", image, CONVERT_MAG);
     
@@ -90,7 +74,7 @@ main (int argc, char **argv)
     
     complex_image_show ("test", image, CONVERT_MAG);
     
-    image_ifft (image);
+    complex_image_ifft (image);
     
     complex_image_show ("test", image, CONVERT_RE);
     
diff --git a/voidcluster.c b/voidcluster.c
new file mode 100644
index 0000000..8ed704a
--- /dev/null
+++ b/voidcluster.c
@@ -0,0 +1,340 @@
+#include <assert.h>
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+#include "image.h"
+
+static int
+X (complex_image_t *image, int i)
+{
+    return (i + image->width) % image->width;
+}
+
+static int
+Y (complex_image_t *image, int j)
+{
+    return (j + image->height) % image->height;
+}
+
+static int
+idx (complex_image_t *image, int x, int y)
+{
+    return Y(image, y) * image->width + X(image, x);
+}
+
+static double
+get (complex_image_t *image, int i, int j)
+{
+    int id = idx (image, i, j);
+
+    return image->red[id].re;
+}
+
+static void
+set (complex_image_t *image, int i, int j, double pixel)
+{
+    int id = idx (image, i, j);
+    image->red[id].re = image->green[id].re = image->blue[id].re = pixel;
+}
+
+static void
+inc (complex_image_t *image, int i, int j, double delta)
+{
+    double v = get (image, i, j);
+    v += delta;
+    set (image, i, j, v);
+}
+
+static void
+find_max (complex_image_t *real, complex_image_t *filtered, int *x, int *y)
+{
+    int i, j, best_x, best_y;
+    double tightest;
+    double f;
+
+    best_x = -1;
+    best_y = -1;
+
+    for (i = 0; i < real->height; ++i)
+    {
+	for (j = 0; j < real->width; ++j)
+	{
+	    if (get (real, j, i) > 0)
+	    {
+		f = get (filtered, j, i);
+		
+		if (best_x == -1 || f > tightest)
+		{
+		    best_x = j;
+		    best_y = i;
+		    tightest = f;
+		}
+	    }
+	}
+    }
+    printf ("cluster: %d %d %f\n", best_x, best_y, tightest);
+
+    *x = best_x;
+    *y = best_y;
+}
+
+static void
+find_min (complex_image_t *real, complex_image_t *image, int *x, int *y)
+{
+    int i, j, best_x, best_y;
+    double loosest;
+    double f;
+
+    best_x = -1;
+    best_y = -1;
+
+    for (i = 0; i < real->height; ++i)
+    {
+	for (j = 0; j < real->width; ++j)
+	{
+	    if (get (real, j, i) == 0)
+	    {
+		f = get (image, j, i);
+		
+		if (best_x == -1 || f < loosest)
+		{
+		    best_x = j;
+		    best_y = i;
+		    loosest = f;
+		}
+	    }
+	}
+    }
+
+    printf ("void: %d %d %f\n", best_x, best_y, loosest);
+
+    *x = best_x;
+    *y = best_y;
+}
+
+static void
+sanity_check (void)
+{
+#if 0
+    int i, j;
+
+    for (i = 0; i < HEIGHT; ++i)
+    {
+	for (j = 0; j < WIDTH; ++j)
+	{
+	    double a = get_filtered (j, i);
+	    double b = compute_filter (j, i);
+
+	    if (fabs (a - b) > 0.0001)
+	    {
+		printf ("warning %d %d got %f computed %f\n", j, i, a, b);
+	    }
+	}
+    }
+
+    printf ("passed sanity check\n");
+#endif
+}
+
+typedef struct
+{
+    int		width;
+    int		height;
+    double	sigma;
+    int		filter_height;
+    int		filter_width;
+    double *	filter_unit;
+} vc_info_t;
+
+static double
+gaussian (vc_info_t *info, int x, int y)
+{
+    double xx = ((x/(double)info->width) * (x/(double)info->width));
+    double yy = ((y/(double)info->height) * (y/(double)info->height));
+
+    return exp (- (xx + yy) / (2 * info->sigma * info->sigma));
+}
+
+static void
+change_filtered (vc_info_t *info, complex_image_t *image, int j, int i, double delta)
+{
+    int k1, k2;
+
+    for (k1 = i - info->filter_height / 2; k1 <= i + info->filter_height / 2; ++k1)
+    {
+	for (k2 = j - info->filter_width / 2; k2 <= j + info->filter_width / 2; ++k2)
+	{
+	    double u = info->filter_unit [(k1 - i + info->filter_height/2) * info->filter_width +
+					  (k2 - j + info->filter_width/2)];
+	    double d = delta * u;
+
+	    inc (image, k2, k1, d);
+	}
+    }
+}
+
+#if 0
+static void
+filter (vc_info_t *info, int j, int i)
+{
+    info.filtered[j][i] = compute_filter (j, i);
+}
+#endif
+
+static double
+compute_filter (vc_info_t *info, complex_image_t *image, int j, int i)
+{
+    int k1, k2;
+    double r;
+
+    r = 0.0;
+    for (k1 = i - info->filter_height / 2; k1 <= i + info->filter_height / 2; ++k1)
+    {
+	for (k2 = j - info->filter_width / 2; k2 <= j + info->filter_width / 2; ++k2)
+	{
+	    double g = gaussian (info, k2 - j, k1 - i);
+
+	    r += g * get (image, k2, k1);
+	}
+    }
+
+#if 0
+    if (get(j, i) == 1)
+	assert (r > 0);
+#endif
+
+    return r;
+}
+
+void
+void_cluster (complex_image_t *image, int n_points)
+{
+    complex_image_t *filtered, *unit;
+    vc_info_t info;
+    int i, j;
+
+    info.width = image->width;
+    info.height = image->height;
+    info.sigma = 1 / sqrt (n_points * M_PI);
+    info.filter_width = 2 * (int)(floor (4 * info.sigma * info.width)) + 1;
+    info.filter_height = 2 * (int)(floor (4 * info.sigma * info.height)) + 1;
+    info.filter_unit = malloc (info.filter_width * info.filter_height * sizeof (double));
+
+    unit = complex_image_new (info.filter_width, info.filter_height);
+
+    double t = 0;
+    /* Initialize unit filter */
+    for (i = 0; i < info.filter_height; ++i)
+    {
+	for (j = 0; j < info.filter_width; ++j)
+	{
+	    double v =
+		gaussian (&info, j - info.filter_width / 2, i - info.filter_height / 2);
+	    
+	    info.filter_unit[i * info.filter_width + j] = v;
+
+	    t += v;
+
+	    set (unit, j, i, info.filter_unit[i * info.filter_width + j]);
+	}
+    }
+
+    printf ("total %f\n", t);
+    
+#if 0
+    complex_image_show ("unit", unit, CONVERT_RE);
+#endif
+
+    
+    /* Make filtered version of image */
+    filtered = complex_image_new (image->width, image->height);
+
+    for (i = 0; i < info.height; ++i)
+    {
+	for (j = 0; j < info.width; ++j)
+	{
+	    if (get (image, j, i))
+		change_filtered (&info, filtered, j, i, get (image, j, i));
+	}
+    }
+
+#if 0
+    complex_image_show ("unfiltered", image, CONVERT_RE);
+
+    complex_image_mul (filtered, 1/128.0);
+    complex_image_show ("filtered", filtered, CONVERT_RE);
+    complex_image_mul (filtered, 128.0);
+#endif
+    
+    /* void/cluster */
+    for (;;)
+    {
+	int cluster_x, cluster_y, void_x, void_y;
+	double delta_cluster, delta_void;
+	double t;
+
+	find_max (image, filtered, &cluster_x, &cluster_y);
+	find_min (image, filtered, &void_x, &void_y);
+
+	delta_cluster = get (image, void_x, void_y) - get (image, cluster_x, cluster_y);
+
+	change_filtered (&info, filtered, cluster_x, cluster_y, delta_cluster);
+
+	delta_void = get (image, cluster_x, cluster_y) - get (image, void_x, void_y);
+	change_filtered (&info, filtered, void_x, void_y, delta_void);
+
+	t = get (image, cluster_x, cluster_y);
+	set (image, cluster_x, cluster_y, get (image, void_x, void_y));
+	set (image, void_x, void_y, t);
+
+    find_max:
+	find_max (image, filtered, &cluster_x, &cluster_y);
+	if ((cluster_x == void_x && cluster_y == void_y) || delta_cluster == delta_void)
+	    break;
+    }
+}
+
+int
+main ()
+{
+    complex_image_t *image;
+    int i, j, k;
+    GdkPixbuf *pb;
+
+#define N_POINTS 4096
+#define W 256
+#define H 256
+
+    g_type_init ();
+
+    srand48 (time (NULL));
+    srand (time (NULL));
+    
+    image = complex_image_new (W, H);
+    for (i = 0; i < N_POINTS; ++i)
+    {
+	do
+	{
+	    j = rand() % W;
+	    k = rand() % H;
+	} while (get (image, j, k) == 1.0);
+
+	set (image, j, k, 1.0);
+    }
+    
+    void_cluster (image, N_POINTS);
+
+    pb = pixbuf_from_complex_image (image, CONVERT_RE);
+
+    pb = gdk_pixbuf_scale_simple (pb, MAX (W, H), MAX (W, H), GDK_INTERP_NEAREST);
+
+    image = complex_image_from_pixbuf (pb);
+
+    complex_image_show ("asdf", image, CONVERT_RE);
+    
+    complex_image_fft (image);
+    
+    complex_image_show ("asdf", image, CONVERT_MAG);
+
+}