#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);

}