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

#include "SDL.h"

#include "lephisto.h"
#include "log.h"
#include "rng.h"
#include "lfile.h"

#include "perlin.h"

#define NOISE_MAX_OCTAVES         128
#define NOISE_DEFAULT_HURST       0.5
#define NOISE_DEFAULT_LACUNARITY  2.

#define LERP(a, b, x)   (a + x * (b - a))
#define CLAMP(a, b, x)  ((x) < (a) ? (a) : ((x) > (b) ? (b) : (x)))

typedef void* noise_t;

/* Used internally. */
typedef struct {
  unsigned char map[256];   /* Randomized map of indexes into buffer. */
  float buffer[256][3];     /* Random 256x3 buffer. */
  /* Fractal stuff. */
  float H;
  float lacunarity;
  float exponent[NOISE_MAX_OCTAVES];
} perlin_data_t;

static perlin_data_t* noise_new(float hurst, float lacunarity);
/* Basic perlin noise. */
static float noise_get(perlin_data_t* pdata, float* f);
/* Fractional brownian motion. */
/* Turbulence. */
static float noise_turbulence(perlin_data_t* noise, float* f, float octaves);
static void noise_delete(perlin_data_t* noise);

static float lattice(perlin_data_t* pdata, int ix, float fx, int iy,
    float fy, int iz, float fz) {

  int nindex;
  float value;

  nindex = 0;
  nindex = pdata->map[(nindex + ix) & 0xFF];
  nindex = pdata->map[(nindex + iy) & 0xFF];
  nindex = pdata->map[(nindex + iz) & 0xFF];

  value = pdata->buffer[nindex][0]  * fx;
  value += pdata->buffer[nindex][1] * fy;
  value += pdata->buffer[nindex][2] * fz;

  return value;
}

#define SWAP(a, b, t)   t = a; a = b; b = t

#define FLOOR(a)        ((int) a - (a < 0 && a != (int)a))
#define CUBIC(a)        (a * a * (3 - 2 * a))

static void normalize(float f[3]) {
  float magnitude;

  magnitude = 1. / sqrtf(f[0]*f[0] + f[1]*f[1] + f[2]*f[2]);
  f[0] *= magnitude;
  f[1] *= magnitude;
  f[2] *= magnitude;
}

static perlin_data_t* noise_new(float hurst, float lacunarity) {
  perlin_data_t* pdata = (perlin_data_t*)calloc(sizeof(perlin_data_t), 1);
  int i, j;
  unsigned char tmp;
  float f = 1;
  for(i = 0; i < 256; i++) {
    pdata->map[i] = (unsigned char)i;
    pdata->buffer[i][0] = RNGF()-0.5;
    pdata->buffer[i][1] = RNGF()-0.5;
    pdata->buffer[i][2] = RNGF()-0.5;
    normalize(pdata->buffer[i]);
  }

  while(--i) {
    j = RNG(0, 255);
    SWAP(pdata->map[i], pdata->map[j], tmp);
  }

  pdata->H = hurst;
  pdata->lacunarity = lacunarity;
  for(i = 0; i < NOISE_MAX_OCTAVES; i++) {
    /*exponent[i] = powf(f, -H);*/
    pdata->exponent[i] = 1. / f;
    f *= lacunarity;
  }
  return (noise_t)pdata;
}

static float noise_get(perlin_data_t* pdata, float *f ) {
  int n[3];     /* Indexes to pass to lattice function. */
  float r[3];   /* Remainders to pass to lattice function. */
  float w[3];   /* Cubic values to pass to interpolation function. */
  float value; 

  n[0] = FLOOR(f[0]);
  n[1] = FLOOR(f[1]);
  n[2] = FLOOR(f[2]);

  r[0] = f[0] - n[0];
  r[1] = f[1] - n[1];
  r[2] = f[2] - n[2];

  w[0] = CUBIC(r[0]);
  w[1] = CUBIC(r[1]);
  w[2] = CUBIC(r[2]);

  /* 
   * This is the big ugly part that is in dire need
   * of optimisation!!!!
   */
  value = LERP(LERP(LERP(lattice(pdata,n[0], r[0], n[1], r[1], n[2], r[2]),
          lattice(pdata,n[0]+1, r[0]-1, n[1], r[1], n[2], r[2]),
          w[0]),
        LERP(lattice(pdata,n[0], r[0], n[1]+1, r[1]-1, n[2], r[2]),
          lattice(pdata,n[0]+1, r[0]-1, n[1]+1, r[1]-1, n[2], r[2]),
          w[0]),
        w[1]),
      LERP(LERP(lattice(pdata,n[0], r[0], n[1], r[1], n[2]+1, r[2]-1),
          lattice(pdata,n[0]+1, r[0]-1, n[1], r[1], n[2]+1, r[2]-1),
          w[0]),
        LERP(lattice(pdata,n[0], r[0], n[1]+1, r[1]-1, n[2]+1, r[2]-1),
          lattice(pdata,n[0]+1, r[0]-1, n[1]+1, r[1]-1, n[2]+1, r[2]-1),
          w[0]),
        w[1]),
      w[2]);
  
  return CLAMP(-0.99999f, 0.99999f, value);
}

static float noise_turbulence(perlin_data_t* noise, float* f, float octaves) {
  float tf[3];
  perlin_data_t* pdata = (perlin_data_t*) noise;
  /* Init locals. */
  float value = 0;
  int i;

  tf[0] = f[0];
  tf[1] = f[1];
  tf[2] = f[2];

  /* Inner loop of spectral construction, where the fractal is built. */
  for(i = 0; i < octaves; i++) {
    value += ABS(noise_get(noise, tf)) * pdata->exponent[i];
    tf[0] *= pdata->lacunarity;
    tf[1] *= pdata->lacunarity;
    tf[2] *= pdata->lacunarity;
  }

  return CLAMP(-0.99999f, 0.99999f, value);
}

void noise_delete(perlin_data_t* noise) {
  free(noise);
}

/* Generate a 3d nebulae map of dimensions w,h,n with ruggedness rig. */
float* noise_genNebulaeMap(const int w, const int h, const int n, float rug) {
  int x, y, z;
  float f[3];
  int octaves;
  float hurst;
  float lacunarity;
  perlin_data_t* noise;
  float* nebulae;;
  float value;
  float zoom;
  float max;
  unsigned int* t, s;

  /* Pretty default values. */
  octaves = 3;
  hurst = NOISE_DEFAULT_HURST;
  lacunarity = NOISE_DEFAULT_LACUNARITY;
  zoom = rug * ((float)h/768.)*((float)w/1024.);

  /* Create noise and data. */
  noise = noise_new(hurst, lacunarity);

  nebulae = malloc(sizeof(float)*w*h*n);
  if(nebulae == NULL) {
    WARN("Out of memory!");
    return NULL;
  }

  /* Some debug information and time setting. */
  s = SDL_GetTicks();
  t = malloc(sizeof(unsigned int)*n);
  DEBUG("Generating Nebulae of size %dx%dx%d", w, h, n);

  /* Start to create the nebulae. */
  max = 0.;
  for(z = 0; z < n; z++) {

    f[2] = zoom * (float)z / (float)n;

    for(y = 0; y < h; y++) {
    
      f[1] = zoom * (float)y / (float)h;

      for(x = 0; x < w; x++) {
 
        f[0] = zoom * (float)x / (float)w;

        value = noise_turbulence(noise, f, octaves);
        if(max < value) max = value;

        nebulae[z*w*h + y*w+x] = value;
      }
    }

    /* More time magic debug. */
    t[z] = SDL_GetTicks();
    DEBUG("   Layer %d/%d generated in %dms", z+1, n,
        (z>0) ? t[z] - t[z-1] : t[z] - s);
  }

  /* Post filtering. */
  value = 1. - max;
  for(z = 0; x < n; z++)
    for(y = 0; y < h; y++)
      for(x = 0; x < w; x++)
        nebulae[z*w*h + y*w + x] += value;

  /* Cleanup. */
  noise_delete(noise);

  /* Results. */
  DEBUG("Nebulae Generated in %dms", SDL_GetTicks() - s);
  return nebulae;
}

/* Generate tiny nebuale puffs */
float* noise_genNebulaePuffMap(const int w, const int h, float rug) {
  int x, y, hw, hh;
  float d;
  float f[3];
  int octaves;
  float hurst;
  float lacunarity;
  perlin_data_t* noise;
  float* nebulae;
  float value;
  float zoom;
  float max;

  /* Pretty default values. */
  octaves = 3;
  hurst = NOISE_DEFAULT_HURST;
  lacunarity = NOISE_DEFAULT_LACUNARITY;
  zoom = rug;

  /* Create noise and data. */
  noise = noise_new(hurst, lacunarity);
  nebulae = malloc(sizeof(float)*w*h);
  if(nebulae == NULL) {
    WARN("Out of memory!");
    return NULL;
  }

  /* Start to create the nebulae. */
  max = 0.;
  f[2] = 0.;
  hw = w/2;
  hh = h/2;
  d = (float)MIN(hw, hh);
  for(y = 0; y < h; y++) {
    f[1] = zoom * (float)y / (float)h;
    for(x = 0; x < w; x++) {
      f[0] = zoom * (float)x / (float)w;

      value = noise_turbulence(noise, f, octaves);

      /* Make value also depend on distance from center. */
      value *= (d - 1. - sqrtf((float)((x-hw)*(x-hw)+(y-hh)*(y-hh))))/d;
      if(value < 0.) value = 0.;

      if(max < value) max = value;

      nebulae[y*w + x] = value;
    }
  }

  /* Post filtering. */
  /*value = 1. - max;
  for(y = 0; y < h; y++)
    for(x = 0; x < w; x++)
      if(nebulae[y*w+x] > 0.)
        nebulae[y*w + x] += value;*/

  /* Clean up. */
  noise_delete(noise);

  /* Results. */
  return nebulae;
}