#include <stdlib.h>
#include <stdint.h>
#include <cairo.h>
#include <gtk/gtk.h>
#include <gdk/gdkpixbuf.h>
#include "fft.h"
#include "image.h"


static const double J1_roots[] =
{
    3.83170597020751,
    7.01558666981561,
    10.1734681350627,
    13.3236919363142,
    16.4706300508776,
};

#define MAIN_LOBE       (J1_roots[0] / M_PI)
#define LOBES           (J1_roots[4] / M_PI)
#define AT_ZERO         (M_PI / 4)

static double
jinc (double d)
{
    if (d == 0.0)
        return AT_ZERO;
    else
        return j1 (M_PI * d) / (2 * d);
}

/* Jinc windowed jinc function */
static double
jinc_jinc (double x, double y)
{
    double d = sqrt (x * x + y * y);
    double rf;

    if (fabs (d) <= LOBES)
        rf = jinc (d) * jinc (d * (MAIN_LOBE / LOBES));
    else
        rf = 0;

    return rf;
}

#define SIZE 1024

static complex_image_t *
make_jinc_image (void)
{
    complex_image_t *result = complex_image_new (SIZE, SIZE);
    int i, j;

    for (j = 0; j < SIZE; ++j)
    {
	for (i = 0; i < SIZE; ++i)
	{
	    complex_t *r = &(result->red[j * SIZE + i]);
	    complex_t *g = &(result->green[j * SIZE + i]);
	    complex_t *b = &(result->blue[j * SIZE + i]);
	    double step = (2 * LOBES) / (double)SIZE;
	    double x = - LOBES + step / 2 + i * step;
	    double y = - LOBES + step / 2 + j * step;
	    double v;

	    v = jinc_jinc (8 * x, 8 * y);

	    v = 0.5 * (v + 1);
	    
	    r->re = g->re = b->re = 0.0;
	    r->re = v;
	    b->re = 0;

	    r->im = g->im = b->im = 0.0;
	}
    }

    return result;
}

static double
sinc (double x)
{
    return sin (M_PI * x) / (M_PI * x);
}

#define LLOBES (3)

static double
lanczos (double x)
{
    if (fabs (x) < LLOBES)
        return sinc (x) * sinc (x / LLOBES);
    else
        return 0.0;
}

static double
lanczos2d (double x, double y)
{
    return lanczos (x) * lanczos (y);
}

static complex_image_t *
make_lanczos_image (void)
{
    complex_image_t *result = complex_image_new (SIZE, SIZE);
    int i, j;

    for (j = 0; j < SIZE; ++j)
    {
	for (i = 0; i < SIZE; ++i)
	{
	    complex_t *r = &(result->red[j * SIZE + i]);
	    complex_t *g = &(result->green[j * SIZE + i]);
	    complex_t *b = &(result->blue[j * SIZE + i]);
	    double step = (2 * LLOBES) / (double)SIZE;
	    double x = - LLOBES + step / 2 + i * step;
	    double y = - LLOBES + step / 2 + j * step;
	    double v;

	    v = lanczos2d (8 * x, 8 * y);

	    v = 0.5 * (v + 1);
	    
	    r->re = g->re = b->re = 0.0;
	    r->re = 1 - v;
            g->re = v;
            b->re = 0;

	    r->im = g->im = b->im = 0.0;
	}
    }

    return result;
}

static complex_image_t *
make_sinc_image (void)
{
    complex_image_t *result = complex_image_new (SIZE, SIZE);
    int i, j;

    for (j = 0; j < SIZE; ++j)
    {
	for (i = 0; i < SIZE; ++i)
	{
	    complex_t *r = &(result->red[j * SIZE + i]);
	    complex_t *g = &(result->green[j * SIZE + i]);
	    complex_t *b = &(result->blue[j * SIZE + i]);
	    double step = (2 * LLOBES) / (double)SIZE;
	    double x = - LLOBES + step / 2 + i * step;
	    double y = - LLOBES + step / 2 + j * step;
	    double v;

	    v = sinc (8 * x) * sinc (8 * y);

	    v = 0.5 * (v + 1);
	    
	    r->re = g->re = b->re = 0.0;
	    r->re = 1 - v;
            g->re = v;
            b->re = 0;

	    r->im = g->im = b->im = 0.0;
	}
    }

    return result;
}

static complex_image_t *
make_jnw_image (void)
{
    complex_image_t *result = complex_image_new (SIZE, SIZE);
    int i, j;

    for (j = 0; j < SIZE; ++j)
    {
	for (i = 0; i < SIZE; ++i)
	{
	    complex_t *r = &(result->red[j * SIZE + i]);
	    complex_t *g = &(result->green[j * SIZE + i]);
	    complex_t *b = &(result->blue[j * SIZE + i]);
	    double step = (2 * LOBES) / (double)SIZE;
	    double x = - LOBES + step / 2 + i * step;
	    double y = - LOBES + step / 2 + j * step;
	    double vj, vs;

	    vj = (jinc (sqrt (8 * x * 8 * x + 8 * y * 8 * y)));
            vs = (sinc (8 * x) * sinc (8 * y));

	    r->re = g->re = b->re = 0.0;
	    r->re = 0.5 * (vj + 1);
	    g->re = 0.0; // 0.5 * (vs + 1);
            b->re = 0.0;

	    r->im = g->im = b->im = 0.0;
	}
    }

    return result;
}

static complex_image_t *
make_box_image (void)
{
    complex_image_t *result = complex_image_new (SIZE, SIZE);
    int i, j;

    for (j = 0; j < SIZE; ++j)
    {
	for (i = 0; i < SIZE; ++i)
	{
	    complex_t *r = &(result->red[j * SIZE + i]);
	    complex_t *g = &(result->green[j * SIZE + i]);
	    complex_t *b = &(result->blue[j * SIZE + i]);
	    double step = (2 * 0.5) / (double)SIZE;
	    double x = - 0.5 + step / 2 + i * step;
	    double y = - 0.5 + step / 2 + j * step;
	    double v;

            if (fabs (8 * x) < 0.5 && fabs (8 * y) < 0.5)
                v = 1;
            else
                v = 0;
                
            
	    r->re = g->re = b->re = 0.0;
	    b->re = v;
	    g->re = 0.0; // 0.5 * (vs + 1);
            r->re = 0.0;

	    r->im = g->im = b->im = 0.0;
	}
    }

    return result;
}


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) / 128.0, 4) + 1)) , arg);
}

static void
low_pass (complex_image_t *image, double d)
{
    int w = image->width;
    int h = image->height;
    int i, j;
    double d2 = d * d;
    int c = h / 2;
    
    for (i = 0; i < h; ++i)
    {
	for (j = 0; j < w; ++j)
	{
	    int idx = i * w + j;
	    double dist = (c - i) * (c - i) + (c - j) * (c -j);
	    double t;
	    
	    image->red[idx] = filter (image->red[idx], dist);
	    image->green[idx] = filter (image->green[idx], dist);
	    image->blue[idx] = filter (image->blue[idx], dist);
	}
    }
}

int
main (int argc, char **argv)
{
    char *input;
    GdkPixbuf *pb;
    complex_image_t *image;
    
    g_type_init ();
    
    image = make_lanczos_image ();
    
    g_print ("width, height %d %d\n", image->width, image->height);
    
    complex_image_show ("test", image, CONVERT_RE);
    
    complex_image_fft (image);
    
    complex_image_show ("test", image, CONVERT_MAG);
    
    complex_image_ifft (image);
    
    complex_image_show ("test", image, CONVERT_RE);
    
    return 0;
}