#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include "fft.h"

#define N_SAMPLES 512

#define K   ((2 * M_PI) / N_SAMPLES)

static double
ssin (double f)
{
    return 0.5 * (sin (K * f) + 1.0);
}

static double
scos (double f)
{
    return 0.5 * (cos (K * f) + 1.0);
}

int
main ()
{
    complex_t *noise = malloc (sizeof (complex_t) * N_SAMPLES * N_SAMPLES);
    int i, j;
    double max, min;

    srand48 (time (0));
    
    for (i = 0; i < N_SAMPLES; ++i)
    {
	for (j = 0; j < N_SAMPLES; ++j)
	{
#if 0
	    noise[i * N_SAMPLES + j].re = 0.6 * (drand48() + drand48()) - 0.2 * drand48();
#endif
	    noise[i * N_SAMPLES + j].im = 0;
	    noise[i * N_SAMPLES + j].re = 0.5 * (drand48() + drand48());
	    noise[i * N_SAMPLES + j].re = 0.6 * (drand48() + drand48()) - 0.2 * drand48();
	}
    }

    fft_2d (noise, N_SAMPLES);
    
    shift_2d (noise, N_SAMPLES);

    max = noise[0].re;
    min = noise[0].re;

    for (i = 0; i < N_SAMPLES; ++i)
    {
	for (j = 0; j < N_SAMPLES; ++j)
	{
	    double d = sqrt ((i - N_SAMPLES/2) * (i - N_SAMPLES/2) +
			     (j - N_SAMPLES/2) * (j - N_SAMPLES/2));
	    double f = (d / (N_SAMPLES)); // * (d / N_SAMPLES);
	    

	    if (d != 0.0)
	    {
		double factor = 0;

#define CHECK(n)							\
		((i <= (n) || i >= N_SAMPLES - ((n) + 1))  ^		\
		 (j <= (n) || j >= N_SAMPLES - ((n) + 1)))
		
#if 0
		if (CHECK(1))
		{
		    factor = 1;
		}
		else if (CHECK(42))
		{
		    factor = 1;
		}
		else
		{
		    factor = 0.1;
		}
#endif
		if (f <= 0.5)
		{
		    factor = pow (0.5 * (cos (M_PI * (1 - 2 * (f))) + 1.0), 8);
		}
		else
		{
		    factor = 1;
		}
#if 0
		noise[i * N_SAMPLES + j].re *= 2.4 * f;
		noise[i * N_SAMPLES + j].im *= 2.4 * f;
#endif
#if 0
		if (f < 0.2)
		    factor = 0;
		else
		    factor = 2.5 * (f - 0.2);
#endif

		noise[i * N_SAMPLES + j].re *= factor;
		noise[i * N_SAMPLES + j].im *= factor; // 8 * N_SAMPLES * N_SAMPLES * factor;
	    }
	}
    }

    printf ("%f %f\n", noise[i * (N_SAMPLES/2) + N_SAMPLES/2].re, noise[i * (N_SAMPLES/2) + N_SAMPLES/2].im);
    
    shift_2d (noise, N_SAMPLES);
    
    ifft_2d (noise, N_SAMPLES);

#define N_BUCKETS  128
    
    int hist[N_BUCKETS];
    int bad;
    bad = 0;
    memset (hist, 0, sizeof hist);
    for (i = 0; i < N_SAMPLES * N_SAMPLES; ++i)
    {
	int v = (noise[i].re * N_BUCKETS);
	if (v == N_BUCKETS)
	    v = N_BUCKETS - 1;
	else if (v >= N_BUCKETS || v < 0)
	{
	    printf ("bad: %f (%d)\n", noise[i].re, v);
	    bad++;
	}
	else
	    hist[v]++;
    }
    int total = 0;
    for (i = 0; i < N_BUCKETS; ++i)
    {
	printf ("%d: %d\n", i, hist[i]);
	total += hist[i];
    }
    printf ("bad: %d (%f %%)\n", bad, 100.0 * bad / (N_SAMPLES * N_SAMPLES));
    printf ("total: %d \n", total + bad);
    
    show_image ("blue noise", noise, N_SAMPLES);
    
    return 0;
}